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INVESTIGATION OF SHOCK-INDUCED 
COMBUSTION PAST BLUNT PROJECTILES 
J. K. Ahuja 1 and S. N. Tiwari 2 

SUMMARY 

A numerical study is conducted to simulate shock-induced combustion in premixed 
hydrogen-air mixtures at various free-stream conditions and parameters. Two-dimensional 
axisymmetric, reacting viscous flow over blunt projectiles is computed to study shock- 
induced combustion at Mach 5.11 and Mach 6.46 in hydrogen-air mixture. A seven- 
species, seven reactions finite rate hydrogen-air chemical reaction mechanism is used 
combined with a finite-difference, shock-fitting method to solve the complete set of 
Navier-Stokes and species conservation equations. In this approach, the bow shock 
represents a boundary of the computational domain and is treated as a discontinuity across 
which Rankine-Hugoniot conditions are applied. All interior details of the flow such as 
compression waves, reaction front, and the wall boundary layer are captured automatically 
in the solution. Since shock-fitting approach reduces the amount of artificial dissipation, 
all the intricate details of the flow are captured much more clearly than has been possible 
with the shock-capturing approach. This has allowed an improved understanding of the 
physics of shock-induced combustion over blunt projectiles and the numerical results can 
now be explained more readily with one-dimensional wave-interaction model than before. 
For Mach 5. 1 1 the flow field is found to be unsteady with regular periodic oscillations 
of the reaction front. There is a progression of higher frequency and lower amplitude 
oscillations as the Mach number is increased with a steady flow observed at some point 
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above the C-J velocity. Numerical results show good qualitative agreement with the 
ballistic range shadowgraphs. In addition, the frequency of oscillations, determined by 
using the Fourier power spectrum is found to be in good agreement with the experiment. 

Various parameters for the triggering of the instabilities have been identified. Projec- 
tile diameter is one of the parameter and an unstable reaction front can be made stable by 
choosing an appropriate small diameter projectile. The other parameter is the heat release 
rate which, in turn, depends upon the free-stream pressure. A number of simulations of 
shock-induced combustion past blunt projectiles in regular and large-disturbance regimes 
are also made at a Mach number of approximately 5 and pressures in the range of 0.1 to 
0.5 atm. For a free-stream pressure of 0.1 atm, the reaction front is steady; at a pressure 
of 0.25 atm, the reaction front develops regular, periodic oscillations. As the pressure 
is increased to 0.5 atm, the oscillations become highly pronounced and irregular. Com- 
bustion with periodic oscillations has been classified as a regular regime and combustion 
with large, irregular oscillations has been classified as a large-disturbance regime. These 
calculations are in agreement with the experimental observations from ballistic-range 
tests. The transition from steady reaction front to regular, periodic oscillations, and then 
to large-disturbance regime is explained by a one-dimensional wave-interaction model. 
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Chapter 1 

INTRODUCTION 

1.1 Importance of Overall Research 

The national commitment to the National Aerospace Plane (NASP) program and 
other hypersonic vehicles such as Trans-Atmospheric Vehicle (TAV) and Aero-assisted 
Orbital Transfer Vehicle (AOTV) have generated renewed interest in hypersonic flows. 
Since these vehicles will rely on air-breathing propulsion, hypersonic propulsion is one 
of the key areas being actively researched. For a successful design of the propulsion 
system to be used for NASP, it is essential to have a clear understanding of the physics 
of mixing and combustion at supersonic speeds in order to develop efficient engines. 
The phenomenon of shock-induced combustion/detonation around bodies travelling at 
hypersonic speeds into combustible mixtures is of great theoretical interest because of 
the need to understand the basic mechanism of combustion instability. These chemical 
instabilities are triggered and sustained via a close coupling between the gas dynamics 
and chemical kinetics and are characterized by a periodic density variation behind the 
shock. It also has a practical application in the development of hypersonic airbreathing 
engines in which supersonic flow is maintained throughout the engine to avoid the losses 
associated with slowing down the fluid. In supersonic combustion ramjet (SCRAMJET) 
engines, strong shocks will most likely occur in local regions of the flow. The presence 
of combustible mixture and shocks create conditions where detonation waves can be 
formed. Thus, it is essential to study the instabilities associated with such conditions to 
have a properly designed engine. The term detonation is applied to the process where a 
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shock and reaction front follow each other very closely and are pressure coupled, while 
shock-induced combustion implies that the shock wave and reaction front are decoupled, 
i.e. the reaction front does not influence the shock directly. In this sense, detonation is 
a limiting case of the shock-induced combustion for the types of flow field considered 
in the present investigations. 

Because of the limited experimental data available in this field, computational fluid 
dynamics (CFD) proves to be a viable tool to better understand the physics of supersonic 
combustion and high-speed flows. The present research will be a good validation tool 
for the numerical codes for high-speed chemically reacting flows. 

1.2 Applications and Motivation 

Conventional ramjets are shown in Figs, la and lb. Figure la shows a supersonic 
ramjet with subsonic combustion. In this conventional ramjet engine, free-stream air 
at high supersonic speeds is compressed to a subsonic Mach number at the entrance 
to the combustor. There are three principal components of the ramjet: (i) the diffuser, 
through which air is admitted to the engine and in which the velocity is reduced and ram 
pressure developed; (ii) the fuel injection system, with which fuel is introduced, vaporized 
and distributed; and (iii) the combustor, which includes a flameholder, the combustion 
zone where heat is released, and a nozzle through which the burned gases are ejected 
rearward at high velocity. Fuel is injected into the combustor, and burning takes place 
in a subsonic stream. It is advantageous over the standard gas turbines in the Mach 
number range of 2 to 5, but is disadvantageous at hypersonic speeds. Figure lb shows a 
hypersonic ramjet with subsonic combustion. In this ramjet free-stream air at hypersonic 
speeds is compressed to a subsonic speed. Slowing from hypersonic to subsonic speeds 
will result in large pressure losses and will cause very high temperature of air entering 
the combustor inlet (much higher than the adiabatic fuel/air flame temperature), resulting 
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Figure 1 Schematics of some ramjet configurations: (a) Supersonic ramjet with 
subsonic combustion; (b) hypersonic ramjet with subsonic combustion; 

(c) hypersonic ramjet with supersonic combustion (Scramjet); 

(d) oblique wave detonation engine where the fuel burns by the 
oblique shocks; and (e) ram accelerator or ramjet in a tube concept. 
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in decomposition of the fuel rather than burning. Therefore, the engine will be a drag 
device rather than a thrust device. 

For an efficient propulsion system at hypersonic speeds, the combustion must take 
place at supersonic speeds, for which two modes of propulsion are being proposed; 
namely, the Scramjet (supersonic combustion ramjet) and Shramjet (shock-induced com- 
bustion ramjet). The Scramjet ([1]— [2]) is an integrated airframe-propulsion concept for 
a hypersonic airplane. A schematic of scramjet is shown in Fig. lc. The entire under 
surface of the vehicle is part of the scramjet engine. Initial compression of the air takes 
place through the bow shock from the nose of the aircraft. Further compression takes 
place inside a series of modules near the rear of the aircraft, thus increasing its pressure 
and temperature. In the combustor, fuel (usually hydrogen) is injected into the hot air 
by a series of parallel and perpendicular injectors where mixing and combustion takes 
place at supersonic speeds. The expansion of burned gases is partially realized through 
nozzles in the engine modules but mainly over the bottom rear surface of the aircraft. At 
high Mach numbers, the fuel and air do not have enough time for mixing and, therefore, 
the combustion efficiency decreases. Thus, in order to get the desired mixing, the length 
of the combustor has to be long. Since the highest pressure and temperature in the en- 
gine occur in the combustor, it has to be very strong; combined with the long length, it 
increases the weight and the drag of the vehicle. 

In order to reduce the size of the combustor, shock-induced combustion (Shram- 
jet [3]) has been proposed, where, a shock is employed to increase the temperature of 
premixed fuel and air to a point where chemical reaction will start. Figure Id shows 
the schematic of Shramjet engine. In Shramjet, fuel is injected well upstream of the 
combustor where temperatures are relatively low and this improves the fuel air mixing. 
Apparent advantages of the Shramjet over the Scramjet engine includes very short-length 
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combustors and simple engine geometries. The Shramjet’s ability to operate at lower 
combustor inlet pressures will allow the vehicle to operate at a lower dynamic pressure 
which lessens the heating loads on the airframe. Shramjet is expected to perform better 
than scramjet in the Mach 12 to 15 range. In another concept which is called the ram 
accelerator or ramjet-in-tube [4], a shaped projectile is fired into a tube filled with a 
premixed gaseous fuel/oxidizer mixture. Figure le shows such a concept. There is no 
propellant on board the projectile. Ignition of the fuel/oxidizer mixture is achieved by 
means of a series of shock waves that increase its temperature until the ignition tempera- 
ture is reached. The resulting energy release develops high pressure behind the projectile 
and this accelerates the projectile to high velocities. The ram accelerator concept has 
the potential for a number of applications, such as hypervelocity impact studies and as 
a mass launcher system. 

1.3. Literature Survey 

In the 1950s the major research effort was directed toward understanding the initia- 
tion, structure, and instability mechanism of the detonation wave. However, most of the 
studies were concerned with the propagation of Chapman- Jouget type detonation in tubes 
filled with combustible mixtures. The main characteristics of these flows is the presence 
of two distinct fronts; the bow shock and the reaction front. These two fronts are sep- 
arated from each other by a distance equal to the induction length. Another interesting 
feature of these flows is the oscillatory behavior of the flow field; the entire flow field 
pulsates in a periodic manner with a characteristic frequency. 

The Chapman-Jouget (CJ) velocity of a mixture (velocity with which a normal det- 
onation propagates in the mixture) is an important parameter and is a characteristic of 
the mixture. If the projectile is travelling at a velocity lower than the CJ velocity, the 
flow field is observed to be highly unsteady and the free-stream velocity is referred to as 
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underdriven, while if the projectile velocity is higher than the CJ velocity of the mixture, 
the flow field appears steady and the free-stream velocity is referred to as overdriven. 

In the past, many researchers have conducted ballistic range experiments to study 
supersonic combustion/detonation. In these experiments, projectiles were fired in dif- 
ferent fuel-air mixtures, and detonation structures around the projectiles were recorded. 
Zeldivich and Shlyapintokh [5] suggested in 1960 that combustion can be stabilized by 
the shock wave produced by bodies moving at supersonic speeds in combustible mix- 
tures. This technique was widely used to study combustion in ballistic range facilities by 
firing projectiles at supersonic speeds into quiescent combustible gas mixtures. In these 
experiments, projectiles were fired in different premixed fuel air mixtures and detonation 
structures around the projectiles were recorded. If the projectile is flying above the C-J 
velocity of the gas mixture, the detonation or reaction front structure shows a coupled 
shock-deflagration system near the stagnation line of the body. These two fronts separate 
from each other as one moves away from the stagnation line. The separation between 
the two fronts occurs as soon as the velocity component normal to the bow shock is 
equal to the detonation velocity. The separation between the bow shock and the reaction 
front is called the induction zone. 

In 1961 Ruegg and Dorsey [6] investigated the problem and effects of stabilizing 
combustion on 20 mm diameter spherical projectiles flying through a quiescent com- 
bustible mixture. Combustion produced detectable effects on the shape and positions of 
shock wave at Mach numbers between 4 and 6.5 and above a pressure of one-tenth at- 
mosphere. Ignition delay was observed behind the bow shock, thus causing a separation 
between the bow shock and the reaction front. Strong combustion driven oscillations 
were also observed with frequencies up to one-tenth megacycles per second. 

Behran et al. [7] conducted similar ballistic experiments by firing 9 mm plastic 
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spheres into hydrogen-air and hydrogen-oxygen mixtures at 1500-3000 ms -1 . They 
also observed that at velocities higher than C-J velocities a steady combustion front is 
established, while at lower velocities unstable forms of oscillations appear. The period 
of oscillations was found to be equal to the induction time for self-ignition. 

Toong and his associates [8-10] conducted a series of experiments using conical and 
spherical projectiles to study the initiation and decay of chemical instabilities. Projectiles 
were fired into lean acetylene- oxygen and stoichiometric hydrogen-air mixtures. They 
proposed the wave interaction model to explain the instabilities in the structure of the 
detonation wave. Their model explains how compression waves can be formed when a 
new reaction front develops in the induction zone between the normal segment of the 
bow shock and the original reaction front. These compression waves lead to a cyclic 
process which is compatible with most of the observed features of the flow. However, the 
strength of the compression waves remained unresolved in their wave-interaction model, 
which is an important factor in determining if such a model is physically possible. Alpert 
and Toong [9] included the effect of the strength of the compression waves and proposed 
a modified form of the wave-interaction model. 

Alpert and Toong [9] investigated detonation-wave structures by firing a spherical 
projectile with a diameter of 12.7 mm at 200 torr initial pressure in a hydrogen-oxygen 
mixture. They proposed that the periodic density variations appear in two main types 
(or regimes) of flow. In the regular regime, the widespread density variations are highly 
regular. The second regime, which is the large-disturbance regime, is characterized both 
by density variations that are less regular but far more pronounced than those of the first 
regime and by distortions of the bow shock. The regular regime has high frequency, low 
amplitude periodic oscillations, while the large-disturbance regime is characterized by 
low frequency, high amplitude oscillations. Alpert and Toong concentrated most of their 
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attention on the large-disturbance regime and proposed a modified form of the wave- 
interaction model to explain the large-disturbance-regime case. 

In 1972 Lehr [11] conducted a detailed experimental study to extend the database for 
a wide range of projectile shapes and combustible mixtures. The projectile shapes tested 
included not only spheres but cones, bi-cones, and flat-nose projectiles. The mixtures 
included hydrogen-air, hydrogen-oxygen, methane-air, and methane-oxygen. The results 
were in general agrement with the previous studies. The shadowgraphs also revealed for 
the first time some three dimensional structure of the flow. 

Oppenheim et al. [12] published an excellent review of detonation research. They 
concentrated on the development of the detonation wave, its stability, and its structure. 
They emphasized that a change to any of the parameters, such as the composition of the 
mixture, its initial pressure, or the diameter of the tube, causes a variation in the amplitude 
of the oscillations of the detonation wave. This variation progressively increases as the 
limits of the detonation are approached in the composition, initial pressure, or diameter. 
Consequently, the wave can possibly become quite unstable. 

Several researchers [13-18] have recently attempted to numerically simulate Lehr’s 
ballistic range experiments [1 1], but have met with limited success. Youngster et al. [13] 
and Lee et al. [14] simulated Lehr’s experimental data for Mach 4.18, 5.11, and 6.46. 
They used Euler equations coupled with species equations to capture the shock and the 
reaction front. The reaction model used was hydrogen-air mixture of six species and an 
inert gas such as Argon or Nitrogen and eight reactions. The flow field was found to be 
steady in contrast to the experimental evidence that the flow field is, indeed, unsteady. 
For the test conditions of stoichiometric hydrogen-air mixture, the detonation wave speed 
of the mixture is Mach 5.11. Experimentally, it has been demonstrated in Lehr’s work 
that Mach 5.11 and 4.18 show structural instabilities of the detonation wave which 
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disappear if the flight Mach number is increased beyond Mach 5.11. Further, the flow 
field was not well resolved. They used 32x32 and 57x41 size grids, respectively, in their 
blunt body calculations. These grids were not sufficient to resolve the flow field correctly. 

Wilson et al. [15] conducted a detailed numerical investigation of the shock-induced 
combustion phenomena. They used Euler equations and a 13-species and 33-reactions 
chemistry model. They showed the validity of the reaction models and the importance 
of grid resolution needed to properly model the flow physics. They did highly resolved 
calculations for Lehr’s Mach 5.11 and Mach 6.46 cases with adaptive grid. The calcula- 
tions were not time accurate, so that the unsteady behavior was not captured. However, 
for cases lower than Mach 5.11, they could successfully capture the instabilities. 

Sussman et al. [16]— [17] also studied the instabilities in the reaction front for a 
Mach number of 4.79. They also used Euler equations and a 13-species and 33-reac- 
tions chemistry model. They have proposed a new formulation based on logarithmic 
transformation. It greatly reduces the number of grid points needed to properly resolve 
the reaction front. They successfully simulated the unsteady case. However, the fre- 
quency was slightly underpredicted. 

Matsuo and Fujiwara [18]— [19] have studied the instabilities of shock-induced com- 
bustion around an axisymmetric blunt body. They used Euler equations and a simplified 
two-step chemistry model. They investigated the growth of periodic instabilities by a 
series of simulations with various tip radii and showed that these periodic instabilities 
are related to shock-standoff distance and induction length. They proposed a new model 
based on McVey and Toong’s model [8]. The instabilities in the reaction front were 
explained by their model. 

The key parameters for the triggering of instabilities has been identified by vari- 
ous parametric studies [20]-[22], Matsuo and Fujiwara [20] and Ahuja and Tiwari [21] 
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showed that an underdriven case, which shows instabilities in the reaction front, can be 
made stable by having an appropriately small size projectile and an overdriven case can 
be made unstable by having a large size projectile. Kumar and Singh [22] concluded that 
the key parameters for triggering these instabilities were projectile velocity, activation 
temperature, projectile nose radius, reaction rate constant, and heat release. 

Tivanov and Rom [23] conducted an analytical study based on an energy equa- 
tion and a chemical rate equation for the flow of a detonable gas mixture over a blunt 
body. They evaluated the conditions for the stability of the detonation process and the 
appearance of the oscillations. The frequency of oscillations matched very well with the 
experimental data. 

Matsuo et al. [24] simulated the regular and large disturbance regime cases of 
Ruegg and Dorsey [6] using two-step chemistry model. With a series of simulations 
the large disturbance regime was explained with a new one-dimensional wave-interaction 
model. Their results revealed that the intensity of heat release was a key parameter in 
determining the regime of the unsteady flow. Flow features of the unsteady combustion 
with low-frequency and high-amplitude oscillation, known as large-disturbance regimes, 
are reproduced when the concentration of the heat release is very high. For moderately 
high heat release, a high-frequency, low-amplitude periodic unsteadiness that belongs to 
regular regimes was observed. 

Ahuja et al. [25]-[29] used the Navier-Stokes equations with a nine-species and 
eighteen-reactions and seven-species and seven-reaction H 2 -air reaction model to simu- 
late Lehr’s Mach 5.11 and 6.46 cases. They used both shock-capturing and shock-fitting 
method to resolve the flow field for the above two cases. Shock-fitting method gave 
much better results with all the intricate flow features very well resolved. The Mach 5.1 1 
case was found to be unsteady while the Mach 6.46 case was macroscopically stable. 
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The frequency of oscillations was found to be in good agreement with the experimen- 
tal frequency. Ahuja et al. [30] also simulated Ruegg and Dorsey’s [6] regular and 
large-disturbance cases and tied their numerical results with Matsuo and Fujiwara’s [24] 
one-dimensional wave-interaction model to explain the large-disturbance regime. 

All the investigations considered so far dealt with ballistic range experiments and 
their simulations. The primary focus was on blunt bodies, i.e., detached detonation waves. 
Attention will now be diverted to attached detonation waves, i.e., flow fields involving 
sharp nosed bodies. The experimental evidence of attached detonation waves is very 
sketchy, and even the experiments themselves are controversial. 

Gross et al. [31] conducted a series of experiments involving shock-induced combus- 
tion/detonation in a H 2 -air mixture. The experiments were conducted in a Mach 3-plus 
flow of H 2 -air mixture at total temperature of 833 K. They claimed to have obtained a 
steady, planer and reproducible attached detonation wave. Pratt et al. [32] re-examined 
their results and concluded that what had been interpreted as an attached oblique detona- 
tion wave was in fact a non-reacting strong oblique shock formed due to shock-induced 
combustion downstream of the test section, off camera, which increased the back pressure 
and, thus, supported the strong shock wave. 

Adelman et al. [33] designed a proof-of principle experiment to be conducted at the 
NASA Ames Research Center. But, they incorrectly identified the C-J turning angle as 
the maximum stabilization angle of the oblique detonation wave. 

In the last few years, a large number of numerical computations involving attached 
detonation wave have appeared in the literature. Fort et al. [34] studied the shock induced 
combustion of premixed H 2 -air over a ramp using the RPLUS code. The chemistry was 
modelled by an 18-step model. The simulations were carried out for various wedge an- 
gles for both viscous and inviscid flows, and it was concluded that any study of attached 
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oblique detonation wave that ignores viscous effects is incomplete, if not in error. Their 
results also showed an ignition hysteresis phenomena. 

Singh et al. [35] conducted a numerical study to address the structural stability of 
an oblique detonation wave. They used Navier-Stokes equations coupled with a seven- 
species and seven-reaction H 2 -air model. The viscous effects were included to account 
for the elliptic influence of the viscosity in the near lip region. All the calculations were 
time accurate. They concluded that an oblique detonation wave is a stable phenomenon 
as long as a sufficient amount of overdrive is present. 

Li et al. [36] studied combustion mechanisms applicable to ram accelerators using 
numerical simulations. Their study showed that it is possible to generate steady deto- 
nation waves over wedge surfaces under appropriate flow and mixture conditions. They 
used the Flux-Corrected Transport (FCT) algorithm to solve mass, momentum, energy 
and species density equations. The chemical reactions were simulated by using a two 
step reaction model. In addition, all diffusive transport processes were neglected. 

The instability in the structure of the reaction front originates in the induction zone 
which separates the bow shock and the exothermic reaction front in the nose region of 
the flow field and then spreads outwards. In order to capture the physical instabilities, 
the calculations must be carried out for long times to ensure that all relevant time scales 
are being captured. Since all numerical schemes have some numerical diffusion, which 
is dependent on the grid resolution, a coarse grid may damp these oscillations. Further, 
the numerical damping added to the scheme in the vicinity of the reaction front may 
damp or alter the instability modes. The objective of this study is to investigate, in 
detail, the shock-induced combustion phenomena for the premixed stoichiometric H 2 -air 
mixture flow at hypersonic speeds. The analysis is carried out using the axisymmetric ver- 
sion of the SPARK2D code [37], which incorporates a 7-species, 7-reactions combustion 
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model for hydrogen-air mixture. The code also incorporates a 9-species and 18-reactions 
hydrogen-air model. The code has both shock-capturing and shock-fitting capabilities. 


13 



Chapter 2 

PHYSICAL PROBLEMS AND CONCEPTS 

In this chapter we shall be explaining some of the shadowgraphs of Lehr [11] and 
Ruegg and Dorsey [6]. Later in the chapter some of the concepts of shock-induced 
combustion shall be discussed followed by various models to explain some of the 
phenomenon in shock-induced combustion. 

To help explain the experimental shadowgraphs, a schematic of shock-induced 
combustion is shown in Fig. 2.1. The figure depicts a supersonic, one-dimensional 
flow of combustible gas which encounters a stationary normal shock. The temperature 
upstream of the shock is too low for combustion but the temperature behind the shock 
is high enough to induce combustion. As the gas passes through the shock, the rise in 
temperature initiates the chemical reactions which leads to burning. Because it takes 
time for the gas to reach the ignition temperature, the gas during that time shall travel 
downstream before ignition can occur if the fluid motion is relatively fast. This cause a 
separation between the bow shock and the energy release front. This separation between 
the bow shock and the energy release front is referred to as the induction zone. The 
induction zone is characterized by an almost constant values for the fluid quantities 
such as temperature, density, and pressure after the shock. The length of the ignition 
zone is determined by the ignition time at the post-shock conditions and the fluid speed 
downstream of the shock. The pressure across the energy release front is nearly constant, 
the temperature rises, and the density drops. After the energy release front, pressure rises, 
while temperature and density relatively remains constant. What we observed in a tube 
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Figure 2.1 Schematic of shock-induced combustion. 


filled with combustible gas mixture can also take place in external flows where a blunt 
projectile is fired in a combustible gas mixture. 

2.1 Experimental Work of H. F. Lehr 

In 1972 Lehr [11] conducted a detailed experimental study for a wide range of 
projectile shapes and combustible mixtures. The experimental shadowgraphs of Lehr 
were chosen to support the numerical simulations presented in this work. This is because 
the quality of the shadowgraph available is excellent and therefore numerical results can 
be more clearly tied to the experimental results. In Lehr’s work spherical nosed projectiles 
of 15 mm diameter were fired in stoichiometric hydrogen-air mixtures for a range of Mach 
numbers so that both steady and unsteady flow phenomenon are represented. In unsteady 
cases only regular regimes are included in his work. Table 2.1 at the end of the chapter 
describes the various free-stream conditions used in Lehr’s work and the corresponding 
frequencies of oscillations. In the current numerical simulations we have chosen two 
cases from the spectrum of cases of Lehr’s work. These are Mach 5.11 and 6.46 cases. 
Reason we choose Mach 5.11, because it is the most controvercial case. Experimentally 
it has been observed that Mach 5.11 shows instabilities in the reaction front whereas 
numerically no researcher has ever been able to show these instabilities. Though for 
lower than Mach 5.11, earliar researchers have successfully simulated the instabilities. 
Mach 6.46 is the superdetonative case where the instabilities disapperas. 

Ballistic range shadowgraph pictures for Mach 5.11 and Mach 6.46 from Lehr’s 
experiments are shown in Figs. 2.2 and 2.3, respectively. In both cases, a free-stream 
temperature of 292 K and a pressure of 42663.2 N/m 2 (320 mm of Hg) is used along 
with a stoichiometric mixture of hydrogen-air. The projectile diameter was 15 mm. At 
these conditions the C-J Mach number of the mixture is 5.11. 

Figure 2.2 shows two discontinuities separated from each other. The outer front 
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Figure 2.3 Shadowgraph of a spherical nose projectile moving at Mach 6.46 into a premixed stoichiometric 

hydrogen-air mixture (Courtesy of Dr. H. F. Lehr). 


is the bow shock and the inner front is the reaction front produced by ignition of the 
heated H 2 -air mixture. The separation between the shock front and the reaction front is 
called an induction zone. The separation between the two is minimum near the stagnation 
point and increases as the shock curves around the body, due to increase in induction 
distance (decrease in post shock temperature) away from the stagnation zone. In general, 
the ignition time is inversely related to the temperature because the chemical reactions 
proceed faster (and the separation becomes smaller) with increasing temperature. A close 
examination of the shadowgraphs reveals that as the flow crosses the bow shock, the color 
changes from light to dark, indicating an increase in density. But, as the flow crosses 
the reaction front, the color changes from dark to light, indicating a decrease in density 
across the reaction front. This is due to a large release of energy across the reaction front, 
causing an increase in the temperature; since the pressure remains relatively constant, the 
density must decrease. Another interesting feature is the presence of corrugation in the 
reaction front. These corrugations are caused by the pulsation of the reaction front. The 
frequency of this pulsation was determined to be 1.96 MHz [11, 38]. 

Figure 2.3 is for the Mach 6.46 case, and it is seen that the reaction front is coupled 
with the shock near the stagnation line and up to about 60-65 degree body angle from 
the stagnation line. This is because of a very high post-shock temperature at Mach 6.46 
that causes the induction zone to become so narrow that it appears that the two fronts are 
merged with each other. Decoupling begins further downstream from the stagnation line 
when the post-shock temperature starts decreasing and, therefore, the induction distance 
increases. Further, both the bow shock and the reaction fronts are smooth without any 
visible instabilities. Thus for an overdriven case of Mach 6.46, the instabilities have 
disappeared. References [11, 38] show other underdriven cases where it has been shown 
that the instabilities in the reaction front become more pronounced as we reduce the 
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projectile velocity lower than the C-J velocity of the mixture. In all these cases the 
projectile diameter was fixed as 15 mm. 

2.2 Experimental Work of Ruegg and Dorsey 

Depending upon the magnitude of various parameters, the detonation wave can be 
stable or unstable. Further, the instabilities in the detonation-wave structure can be highly 
periodic; these instabilities are termed as “regular regime”. The instabilities can also be 
highly pronounced and irregular with large-amplitude oscillations and a distorted bow 
shock. These oscillations have been classified in the literature as the “large-disturbance 
regime”. Ruegg and Dorsey [6] fired a projectile with a diameter of 20 mm in a 
stoichiometric hydrogen-air mixture at various free-stream pressures and Mach numbers. 
Tests were made at pressures from 0.026 to 1 atm (absolute) and at Mach numbers up to 
6.5. The effects of combustion on the wave shape and position were detected at pressures 
of 0.1 atm and higher and for a range of Mach numbers from 4 to 6.5. Results between 
pressures of 0.5 and 1 atm were qualitatively similar. When the Mach number was kept 
approximately constant in the range of 4.8 to 5 and only the free-stream pressure was 
changed from 0.1 to 0.25 atm, the reaction front changed from a stable to a periodic 
unstable front. For the same Mach-number range, the reaction front instabilities became 
highly pronounced (with large amplitudes) and irregular when the pressure was increased 
to 0.5 atm. The bow shock was completely distorted by the reaction front, and the period 
of oscillation became much higher than the regular periodic case. 

In the present study, some of the cases of Ruegg and Dorsey’s [6] experimental 
work have been simulated numerically to demonstrate the transition from steady regime 
to regular unsteady regime to large-disturbance unsteady regime as the free-stream 
pressure is changed from 0.1 to 0.25 to 0.5 atm (the free-stream Mach number is kept 
approximately constant around 5). Table 2.2 at the end of this chapter describes the 
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Figure 2.4 Ruegg and Dorsey’s ballistic experiments shadowgraphs (a) M=5.4; p=0.5 atm; Air (b) M=43; 
p=0.5 atm, (c) M=4.8; p=0.5 atm (d) M=5.5; p=0.5 atm (e) M=6.3; p=0.5 atm and (f) M=5.0; p=0.1 atm. 



Figure 2.5 Ruegg and Dorsey’s ballistic experiments shadowgraphs for p=0.25 atm 
(a) M=4.4; .3H 2 + .7N 2 (b) M=4.5 (c) M=4.9 (d) M=5.1 (e) M=5.9 and (f) M=6.5. 



Figure 2.6 Enlarged experimental shadowgraph for M=5.0; p=0.1 atm. 



Figure 2.7 Enlarged experimental shadowgraph for M=4.9; p=0.25 atm. 




Figure 2.8 Enlarged experimental shadowgraph for M=4.8; p=0.5 atm 


various free-stream conditions used in Ruegg and Dorsey’s work. Figures 2.4 and 2.5 
show schlieren shadowgraphs of a series of experiments done by Ruegg and Dorsey 
[6] under various free-stream Mach-number and pressure conditions. Figure 2.4a is the 
schlieren shadowgraph of the missile in air; Fig. 2.5a is the shadowgraph for the same 
missile fired in a mixture of nitrogen and hydrogen. These results are presented here 
for comparison with those in the combustible gas. Figures 2.4 and 2.5 show additional 
results for the stoichiometric mixture of H 2 and air. At M = 4.3 and p = 0.5 atm. Fig. 
2.4b shows a discrete region of reaction front that surrounds the sphere. A smooth bow 
shock is also visible. As the Mach number is increased, the reaction front and shock 
front merge to give smooth and continuous discontinuities at M = 6.3, as shown in Fig. 
2.4e. However, at a pressure of 0.5 atm and M = 4.8 as shown in Fig. 2.4c and at M = 
5.5 as shown in Fig. 2.4d, the bow shock and the reaction fronts are completely distorted. 
The oscillations in the reaction front become pronounced due to high heat-release rates. 
Figure 2.4f shows the results for M = 5.0 and p = 0.1 atm. At this low pressure, the 
bow shock and the reaction front are both smooth and are separated from each other by 
a large induction region. 

When p = 0.25 atm (Figs. 2.5b-2.5f), a dark combustion wave is clearly visible. 
For M = 4.5 and p = 0.25 atm, Fig. 2.5b clearly shows the smooth bow shock and 
the reaction front separated near the stagnation region; this separation increases as we 
move away from the stagnation region. The key point to' emphasize here is that both the 
bow shock and the reaction front are smooth, with no visible oscillations. An increase 
in Mach number to 4.9, as shown in Fig. 2.5c, triggers periodic instabilities. The bow 
shock is nearly smooth; however, the reaction front appears highly periodic. A further 
increase in Mach number from M = 4.9 (Fig. 2.5c) to M = 5.1 (Fig. 2.5d) to M = 6.5 
(Fig. 2.5f), causes these instabilities to disappear. Consequently, a smooth reaction front 
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evolves, which merges with the bow shock. The separation between the bow shock and 
the reaction front (i.e., the induction zone) decreases continuously in Figs. 2.5d through 
2.5f as the Mach number increases. This decrease results from the high post-shock 
temperature, which causes the two fronts to nearly merge. 

In the present study, we simulate three cases. The first case corresponds to Fig. 
2.4f at Mach number of 5 and pressure of 0.1 atm. This case is steady; the reaction 
and shock fronts are smooth. The second case corresponds to Fig. 2.4c at M = 4.8 and 
p = 0.5 atm. This case is a large-disturbance regime case; it shows a distorted bow 
shock and large-amplitude nonperiodic oscillations of the reaction front. The last case 
is shown in Fig. 2.5c; here, p = 0.25 atm and M = 4.9. In this case, the bow shock is 
separated from the reaction front by an induction zone. The bow shock is quite smooth, 
but the reaction front shows periodic oscillations. This case is a regular- regime case. 
Figures 2.4f, 2.5c, and 2.4c are enlarged for clarity in Figs. 2.6, 2.7 and 2.8 respectively. 
Wave-detachment distances for the shock wave and the combustion wave are compared 
with experimental and analytical results. The results are interpreted with the concepts 
developed by Oppenheim et al. [12] and the model developed by Matsuo et al. [24], 
Simulations are carried out with the shock-fitting technique. 

23 McVey and Toong’s One-Dimensional Wave-Interaction Model 

In order to explain the instabilities in the reaction front both in Lehr’s Mach 5.11 
case and Ruegg and Dorsey’s regular regime case, McVey and Toong proposed a one- 
dimensional wave interaction model. An x-t diagram of the complete cycle of events for 
the postulated instability model is shown in Fig. 2.9. The figure contains the features along 
the stagnation streamline in time. The steps referred to in the following are designated in 
the figure. The cycle of events starts at a time when the contact discontinuity approaches 
the original reaction front. At Step 2 the hot gases behind the contact discontinuity begins 
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Figure 2.9 McVey and Toong’s one-dimensional wave-interaction model. 



to react, generating compression waves which propagates upstream and downstream. At a 
somewhat later time (Step 3), the contact discontinuity reaches the position of the original 
reaction front, extinguishing the reaction at this point and generating rarefaction waves. 
The reaction front begins to recede because of the increasing induction time of the colder 
fluid within the entropy zone. At a later time (Step 4), the compression wave generated 
earlier at the new reaction front interacts with the bow shock, thus strengthening it, 
reflecting a weak rarefaction, and producing another contact discontinuity. The incident 
rarefaction generated by the extinguishing of the original reaction front penetrates the 
bow shock, thus weakening it, and generating a zone of decreasing entropy. The cycle 
of events is completed as the contact discontinuity followed by the zone of decreasing 
entropy approaches the receded reaction front (Step 1). 

McVey and Toong described how the compression waves and contact discontinuities 
present in the interaction model can explain the various features observed in experimental 
shadowgraphs. 

The interaction model of the Alpert and toong added more complexity to the 
model of McVey and Toong by accounting for the strength of the compression waves 
created by the exothermic reactions in the new energy-release front. Accounting for such 
effects was deemed necessary to investigate the large-disturbance regime because it had 
been observed that the oscillations in this regime were more likely to occur when the 
Damkohlar parameter (the ratio of chemical energy released to the sensible energy present 
before reaction) was large. The Alpert and Toong mechanism theorizes that four periods 
of a type similar to the model by McVey and Toong occur within each large-disturbance 
period. Each of the McVey and Toong type of interactions are slightly different and occur 
in such a way that one of them is re-enforced and effects the flow much more than the 
others. This could cause the irregular behavior observed in the large-disturbance case. 
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Because of the complexity of the Alpert and Toong mechanism we shall not discuss this 
method here. Instead we shall be using another simple model developed by Matsuo and 
Fujiwara [24] to explain the large-disturbance regime case. 

2.4 Matsuo and Fujiwara’s One-Dimensional Wave-Interaction Model 
The wave interaction in the regular regime case of Lehr or Ruegg and Dorsey showed 
a completely different structure from that of the large-disturbance regime of Ruegg and 
Dorsey. The regular corrugated structure of regular regime disappears in the large- 
disturbance regime case and is replaced by a non-periodic oscillations of large amplitude. 
The oscillations become much more pronounced with low frequency and high amplitude. 
Further the period of oscillations for large-disturbance regime is 4-5 times that of the 
regular regime case. The mechanism of the large disturbance regime is also dominated 
by the wave interactions but the role of wave and interactions change from those of the 
regular regime. 

The most important point of the mechanism of the regular regime is different 
induction time before and behind the contact discontinuity. On the other hand, the 
large-disturbance regime shows a new feature of periodicity. The extremely strong 
exothermicity occurs on the reaction front, and causes the strong reaction shock toward 
the bow shock and the body surface. The reaction shock is so strong that the gas behind 
the reaction shock is compressed very much, and the exothermic reaction follows and 
accelerates the reaction front. The phenomenon is considered to be onset of “ explosion 
within an explosion”, producing two strong shock waves in opposite directions. The 
forward shock is referred to as “superdetonation” and moves into unbumed gas. In the 
opposite direction a shock moves into the burned gas and is known as “retonation”. 
The mechanism is usually observed as a typical example of deflagration-to-detonation 
transition (DDT). 
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Figure 2.10 Matsuo and Fujiwara’s one-dimensional wave-interaction model 


The detonation wave in the unbumed gas behind the bow shock overtakes and 
penetrates the bow shock, then the reflected shock and the contact discontinuity are 
generated. The intersection point of the bow shock, the detonation wave and the reflected 
shock is called the triple point which is usually observed in the detonation wave structure. 
After the penetration, the detonation wave cannot develop the self sustained detonation 
in front of the projectile body. Eventually the bow shock wave accelerated by the 
penetration of the detonation wave with respect to the projectile body is decelerated, 
and the transition from detonation to shock-deflagration system, which is the ordinary 
shock-induced combustion appears. 

Figure 2.10 is the model proposed by Matsuo and Fujiwara to explain the large- 
disturbance regime. The steps to be referred in the following are indicated by the 
bracketed numbers in the left hand margins of Fig. 2.10. The beginning of the cycle 
is shown at Step 1, when “explosion within an explosion” occurs on the reaction front, 
and then the reaction shocks propagate upstream and downstream. The forward shock 
referred to as “superdetonation” moves into the unbumed gas, and the backward shock 
referred to “retonation” moves into the burned gas. The superdetonation speed is much 
faster than the retonation speed due to acceleration by the following reaction front. At 
Step 2, the detonation wave overtakes and penetrates the bow shock. Then the bow 
shock and the detonation wave creates a triple point and generates a reflected shock and 
a contact discontinuity. The bow shock is accelerated by the penetration, and the gas 
behind the bow shock is much compressed. At Step 3, the bow shock is decelerated and 
the bow shock and the reaction front becomes separated. At Step 4, the retonation wave 
and the reflected shock reach the body surface and the reflected compression waves go to 
the bow shock. The reflected compression wave interacts with the bow shock at Step 5, 
and the contact discontinuity is created. The temperature behind the contact discontinuity 
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is higher than that before it, so that the induction time becomes shorter. At Step 6, the 
contact discontinuity reaches the original reaction front, and the “explosion within an 
explosion” occurs on it and the cycle of events is completed. 
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Table 2.1 H. F. Lehr’s Experimental Ballistic-Range Data 

Projectile diameter =15 mm 
Free-stream pressure = 0.42 atm. (320 mm of Hg) 
Free-stream temperature = 292° K 
Combustible gas: stoichiometric hydrogen-air mixture 
2H 2 + 0 2 + 3.76N 2 -*• 2H 2 0 + 3.76N 2 


Mach No. 

Free-stream velocity, 
u (m/sec) 

Free-stream velocity (u) / m \ 
Detonation velocity (D) Vsec / 

Frequency 

f(MHz) 

6.46 

2605 

1.27 

— 

5.11 

2058 

1.00 

1.96 

5.04 

2029 

0.99 

1.04 

4.79 

1931 

0.94 


4.18 

1685 

0.82 

0.15 | 


Table 2.2 Free-Stream Conditions for Ruegg and Dorsey’s Data 

Projectile diameter = 20 mm 
Combustible gas: stoichiometric hydrogen-air mixture 
2H 2 + 0 2 + 3.76N 2 -* 2H 2 0 + 3.76N 2 


Mach No. 

Poo, atm 

n 

Case 

5.0 

0.1 

300 

Steady Regime 

4.9 

0.25 

300 

Regular Unsteady Regime 

4.8 

0.5 

300 

Large-Disturbance Unsteady 
Regime 
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Chapter 3 

GOVERNING EQUATIONS AND METHOD OF SOLUTION 

In this chapter we shall be describing the basic governing equations used to numerical 
simulate the physical problem for the chemically reacting flows. Also detailed chemistry 
and thermodynamics models used shall be discussed. Finally, to include the effects of 
diffusion of momentum, energy and mass, kinetic theory based diffusion transport models 
are incorporated into the program. Details of the models shall be discussed. 

3.1. Basic Governing Equations 


The physical model for analyzing the flow field is described by the Navier-Stokes and 
species continuity equations. For two-dimensional planar or axisymmetric flows, these 


equations are expressed in physical coordinates as 

<9U dF dG 
dt + dx + dy 


( 3 . 1 ) 


where vectors U, F, G, and H are written as 



pu 

pu 2 - <T X 
pUV - T xy 

F = ( pE - a x )u - T X yV + q x 

pfi(u + Ui ) 
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G = 


pv 

PUV - T X y 

pv 2 - (Ty 

(P-E T X yU (]y 

pf t {v + Vi) 


For axisymmetric flow, 



pv 

(pVU + T X y) 

PV 2 + Tyy - TQ6 

(, pE + p + Tyy ) V + T xy U + q y 

OJ, 


and for planar flow 



'0 ‘ 
0 
0 
0 
uii 


The other terms that appear in the vectors F, G, and H are defined as 


(7 x 


du 

-p + 2p— + AV.u 

ox 


(3.2) 


° y 


dv 

-p + 2 p— + A V.u 


(3.3) 


T xy 


( du dv 

^\dy dx 


) 


(3.4) 


T yy ~ 


2 ( dv v 

3 \ dy y 



(3.5) 
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2 ( v dv du 

T ee = — 3 ^ y y ~ ~dy~ ~dx 


Ul i r ~ 

k Jii + P2^ h 'f' u ' 


-A:— + /g 


pw'Lw. 


2 I „,2 


v-> , „ P u -\- v 

£ = E*./.-! + ^r 


(3.10) 


In Eq. (3.1), only (A^-l) species equations must be considered in the formulation 
because the mass fraction of the species is prescribed by satisfying the constraint equation 


5> = ‘ 


(3.11) 


The binary diffusion equation for the diffusion velocity of the i th species 


u, = u,i + u 3 j 


(3.12) 


vx < = E (^)<“> - “'■) + (/■ - *)(y 


(3.13) 


Note that this equation must be applied only to (N s -l) species. The diffusion velocity 

N, 

for the remaining species is prescribed by satisfying the constraint equation^ f t u t = 0, 

t=i 

which ensures consistency. In Eq. (3.13), we assume that the body-force vector per unit 


37 



mass is negligible. In addition, thermal diffusion is assumed considered to be negligible 
in comparison with the binary diffusion coefficient. 

3.2. Thermodynamic Model 

In order to calculate the thermodynamic quantities, the specific heat for each species 
is first calculated by a fourth-order polynomial in temperature: 

% = A + BiT + C t T 2 + D t T 3 + E X T 4 (3. 14) 

tii 

The coefficients A u B u Ci, D t , and E t for each species are found by a curve fit of the 
data tabulated in Ref. 39. Once we know the specific heat of each species, the enthalpy 

of each species can then be found from 

T 

hi = h? + Jc p ,dT (3.15) 

J'R 

and then the total internal energy can be calculated from Eq. (3.10). In order to determine 
the equilibrium constant (which we shall be needing in the next section) for each chemical 
reaction being considered, we need the information about the Gibbs energy of each 
species. For a constant pressure process, C p /R from the above Eq. (3.14) is first integrated 
over temperature to define the entropy of the species, and the resulting expression 
is integrated again over temperature to obtain the following fifth-order polynomial in 
temperature for the Gibbs energy of each species: 

f = )t‘-(|)t» + F,-G,T (, = 1, JV.) 

(3.16) 

Coefficients F, and G, are defined in Ref. 39. The Gibbs energy of reaction can then 
be calculated as the difference between the Gibbs energy of product species and reactant 
species. 

N a N , 

AG R, = E v ]' 9i - v 'j' 9i (j = l ' -• iVr ) ( 3 - 17) 

•=i i=i 
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The equilibrium constant for each reaction is then found from 


1 ^ f-AG Rj 

w) ex H^r- 


0 = 1 ,-JVr) 


where 

n. n 3 

An j = “ 1l u 'i' = X ’ - Nr ) ( 

1 = 1 1 = 1 

and is the change in the number of moles when going from reactants to products. 


3.3 Chemistry Model 


In the current simulations, a finite-rate chemical reaction of gaseous hydrogen fuel 
and air has been used. That reaction is modeled by two different ways. The reactions are 
first modeled by 9-species, 18-reaction. Next the reactions are modeled by 7-species 
and 7-reactions. Both the models are shortened versions of the Jachimowski’s original 
13-species and 33 reactions hydrogen-air model. In the nine species model eight of the 
chemical species (H 2 , 0 2 , H 2 0, OH, H, O, H0 2 , H 2 0 2 ) are active, and the ninth (N 2 ) 
is assumed inert. 


Chemical-reaction-rate expressions are usually determined by summing the contri- 
butions from each relevant reaction path to obtain the total rate of change of each species. 
Each path is governed by a law-of-mass-action expression in which the rate constants 
can be determined from a temperature-dependent Arrhenius expression. The forward rate 
for each reaction j is determined from the modified Arrhenius law 

*/)= Aj r “' ex p(^j) O' = (3.20) 

The appropriate constants A j,aj, and e } for the H 2 -air reaction system can be found in 
Table 4.1. Knowing the forward rate, and using the equilibrium constant determined in 
the previous section, the backward rate can be defined by 


K bj = ^ = 

K eq 3 


(3.21) 
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Once we have determined the forward and reverse reaction rates, the production rates of 
the eight species (for 9-species and eighteen reaction model) or 6 species (for 7-species 
and 7-reaction model) can be found from the law of mass action. For the general vV r -step 
chemical reaction 


N, N . 

5>;,S,^5>£S,0' = l,..J v r ) (3.22) 


t= 1 t=l 

the law of mass action states that the rate of change of concentration of species i by 


reaction j is given by 


N a 


N, 


Ci = wii - td * fi n n c ™ m 


(3.23) 


777 = 1 


m=l 


All third-body efficiencies are assumed to be equal to 1.0. The net rate of change in 
concentration of species i by reaction j is then found by summing the contributions from 
each reaction. 


Ci = 



(3.24) 


Finally, in the vector H, in Eq. (3.1) the term represents the net rate of production 
of species i in all chemical reactions and is given by 

ui = CM (3.25) 

3.4 Diffusion Models 

Models for the coefficients governing the diffusion of momentum, energy, and mass 
are now determined. The individual species viscosities are computed from Sutherland’s 
law. 


jj_ _ f T \ 3 / 2 T q + S 

Vo \To/ 


40 


T + S 


(3.26) 



where fi o and To are reference values and S is the Sutherland constant. All three values 
are tabulated for the species in Refs. 40 and 41. Once the viscosity of each species has 
been determined, the mixture viscosity is determined by Wilke’s law (Ref. 42) 


Mi 

'‘“ = 2 . T. 

1=1 1 + £ £ 


where 


(72) I 1 + (f )] ' 

The species thermal conductivities are also computed from Sutherland’s law 


(3.28) 


T \ T' 0 + S' 

T'J T + S' 


(3.29) 


but with different values of the reference values ko and T 0 , and the Sutherland’s constant 
1 

S . These values are also taken from Refs. 40 and 41. The mixture thermal conductivity 
is computed by using conductivity values for the individual species and Wassiljewa’s 
formula (Ref. 43), 


iV, . 

“■» = £ x — 

,_1 1 + T, E 


(3.30) 


where - l.O650,y and (fri-j is taken from Eq. 3.28. 

For dilute gases, Chapman and Cowling used kinetic theory to drive the following 
expression for the binary diffusion coefficient D y between species i and j (Ref. 40): 


0.001858T 3 / 2 [^jrjr 


Dy = 


p<7 2 .ft D 


Here, the diffusion collision integral flo is approximated by 


flo = T 


* — 0.145 


+ (T* +0.5)- 
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where 


T* 



Values of the effective temperature T e and effective collision diameter a are taken to be 
averages of the separate molecular properties of each species, giving (Ref. 40) 


+ °i) 

and 

T<„ = (T,,T r ,) 1/2 

Once the binary diffusion coefficients for all species combinations are known, the 
diffusion velocities of each species can be computed from Eq. (3.13). The diffusion 
velocity of each species is the species velocity due to all diffusion processes algebraically 
added to the convection velocity. When computing the diffusion velocities, it is assumed 
as suggested in Ref. 44, that the thermal diffusion coefficient Dr is negligible compared 
with the binary diffusion coefficient. The solution of Eq. (3.13) requires solving a 
simultaneous equation system, with the number of equations equivalent to the number of 
species present for each component of diffusion velocity. 
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Chapter 4 

THE COMBUSTION PROCESS OF HYDROGEN-AIR SYSTEM 

The combustion process involves the oxidation of constituents in the fuel that are 
capable of being oxidized, and can therefore be represented by a chemical equation. 
During a combustion process the mass of each element remains the same. Thus, writing 
chemical equation and solving problems involving quantities of the various constituents 
basically involves the conservation of mass of each element. A brief review of this 
subject, particularly as it applies to the combustion of hydrogen-air system is presented 
in this chapter. 

4.1 Hydrogen-Air Reaction Mechanism 

In hydrogen-air combustion, hydrogen is oxidized by oxygen and chemical energy 
is released and water is formed as product of the reaction. It should be pointed out 
that in the combustion process there are many intermediates products formed during the 
chemical reaction. An elementary chemical reaction is one that takes place in a single 
step. For example, a dissociation reaction such as 

0 2 + M -»• 20 + M (4.1) 

is an elementary reaction because it literally takes place by a collision of an 0 2 molecule 
with another collision partner, yielding directly two oxygen atoms. On the other hand, 
the reaction 


2H 2 + 0 2 2H 2 0 


(4.2) 
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is not an elementary reaction. Two hydrogen molecules don’t come together with one 
oxygen molecule to directly yield two water molecules. Instead Eq. (4.2) is a statement 
of an overall reaction that actually takes place through a series of elementary steps. 


h 2 

2H 

(4.3a) 

0 2 - 

20 

(4.3b) 

h + o 2 - 

OH + 0 

(4.3c) 

o + h 2 - 

OH + H 

(4.3d) 

OH + H 2 -+ 

h 2 o + h 

(4.3e) 


Equations (4.3 a-e) constitutes the reaction mechanism for the overall reaction given by 
Eq. (4.2). Each of Eqs. (4.3a-e) is an elementary reaction. 

The assumption that air is 21.0 per cent oxygen and 79.0 percent nitrogen by volume 
leads to the conclusions that for each mole of oxygen, 79.0/21.0 = 3.76 moles of nitrogen 
are involved. Therefore, when the oxygen for the combustion of hydrogen is supplied 
as air, the overall reaction can be written as 

2H 2 + 0 2 + 3.76N 2 -*• 2H 2 0 + 3.76N 2 (4.4) 

The minimum amount of air that supplies sufficient oxygen for the complete combustion 
of all the hydrogen, is called the “theoretical air” or “stoichiometric air” and the 
combustion mixture is called “stoichiometric hydrogen-air mixture”. The nitrogen acts 
as an inert gas or a diluent. 

There is a broad category of reaction mechanisms that involves homolysis of covalent 
bonds with the production of intermediates posessing unpaired electron called radicals 
(or free radicals). 

A:B — ► A. + B. 
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Almost all small radicals are short-lived, highly reactive species. When they collide 
with other moleculles they tend to react in a way that leads to pairing of their unpaired 
electrons. One way they can do this is by extracting an atom from another molecule. 
For example in the reaction N 2 + O — ► NO + N, the nitrogen radical extracts an oxygen 
atom from oxygen molecule and gives another radical of oxygen. This behavior is 
charecteristic of radical reaction and consists of three basic steps. The first step is called 
the chain-initiating step. In this step radicals are created. The second step is called 
chain-propagating step. In chain-propagating step one radical generates another. The 
third step is the chain-terminating step. This last step occurs less frequently but occurs 
often enough to use up one or both of the reactive intermediates. These intermediate 
process can have time scales similar to the fluid time scales, therefore, fluid dynamic 
simulations requiring finite rate chemistry must account for this. Further it is important 
that all pertinent reactions that may affect the rate process must be included. Values 
of the rate constants for high-temperature hydrogen-air mixture are readily available in 
the literature. But there is always some uncertainty in the published rate constants; 
they are difficult to measure experimentally, and very difficult to calculate accurately. 
Thus the ideal way to create a reaction mechanism is to assemble the important species 
in the hydrogen-air reaction. The most important species in the hydrogen-air reaction 
mechanism are the eight species ( N 2 , O 2 , H 2 , OH, HO 2 , H 2 O, O, H). These species 
form the core of hydrogen air combustion mechanism. 

4.2 Hydrogen-Air Reaction Model 

The hydrogen-air combustion mechanism used in this work is based on the Jachi- 
mowski hydrogen-air model [45], which use 13-species and 33-reactions and is given in 
Table 4.1. The values of reaction- rate constant, temperature coefficient in reaction rate 
expression and activation energy for the various reactions are also tabulated. The current 
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numerical work uses a modified form of this model by using only the most important 
species. We have used both 9-species 18-reactions and 7-species 7-reactions models, 
both of which are tabulated in Tables 4.3 and 4.4 respectively. There is no significant dif- 
ference in the results by these two models. Thus it saves considerably in computational 
time by using the shortened model. 

4.3 Third Body Reactions 

When simple radicals recombine to form a product the energy liberated in the process 
is sufficiently great to cause the product to decompose into the original radicals. 

A + B?= i C + D + energy 

This energy must be removed by a third body M in order to stop recombination. If 
the molecule formed (like C and D above) in a recombination process has a large number 
of internal (generally vibrational) degrees of freedom, it can redistribute the energy of 
formation among these degrees and a third body is not necessary. 

In the above Table 4.3 for 18 reactions model, reactions (6), (7), (8) and (18) needs 
a third body. In some cases the recombination process can be stabilized by the formed 
molecules radiately dissipating some energy or colliding with a surface and dissipating 
energy in this manner. 

Now let us consider the third body reaction 

* 7 , 

A + B + M, *5 C + M, (4.4) 

k k 

where Mi can be any species present in the fluid. Here M,- plays the role of a catalyst 
and as a result the forward and backward reaction rates will depend on which species 
is involved as the third body. 

Rate of formation of species C is 

= K f> C A C B C Ui (4.5) 
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Rate of consumption of species C is 


dCg 

dt 


Kb,CcCu, 


Therefore, net rate of formation of species C is 

- K/.CaCbCm, - K bt C c C M , 


At equilibrium we have 


and Eq. (4.7) becomes 


dC c 

dt 


= 0 


KMCaUCbUCmX = I<b.(Cc) e (C Mt ) e 


(4.6) 


(4.7) 


(4.8) 


or 


IC = 


I< 


/. 


( C c) e 


I<b, (CaUCb)' 

where subscript e stands for equilibrium and K e is the equilibrium constant. Note that 
this is independent of the catalyst species M,-. We can rewrite Eq. (4.7) as 


dC c 
dt 


Kf,C u , 


CaCb - jrC c 


(4.9) 


Now since M, could be any species, so Eq. (4.9) actually represents N reactions (one for 
each of the N species present in the fluid). Total production of species C is the sum of 
the production of C from all of the reactions. 

Therefore, 


dC c 

dt 

dCc _ 
dt 


N r / i xi 

= E[% C m.(CaCb--^C c J 

N 

E( a '/, c m ,) 

L i=i 


C A C B - -^rCc 


(4.10) 
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4.4 Third Body Efficiencies 

We can relate the reaction rate co-efficients, Kf, corresponding to the various species 
acting as third bodies to the reaction rate of one of the species acting as a third body. 
For example, if we choose nitrogen as the third body we can write 

* = (*/)* 

where r/, is the third body efficiency of the ith species. In the above equation i could be 
any species and in the denominator instead of nitrogen it could be any particular species 
as reference. 

Therefore, we can write Eq. (4.10) as 

^ = (*/)„, (c a Cb - ^Cc) < 411 > 

Table 4.2 shows some third body efficiencies for several reactions where H 2 and H 2 0 
are the third bodies, i.e., the collision partner denoted by M in some of the chemical 
equations. 
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Table 4.1 Jachimowski’s Hydrogen-Air Model 


EZ 

Reaction 

rz — 

a j 

e j 


H 2 + 0 2 = OH + OH 


0. 



OH + H 2 =H 2 0 + H 


0. 

Hi 


H + 0 2 = OH +0 

2.20E14 

0. 



0 + H 2 = OH + H 

1.80E10 

1. 



OH + OH = H 2 0 + 0 

6.3E12 

0. 

1090. 


H + OH = H 2 0 + M 

2.20E22 

-2. 

0. 

7 

H + 0 = OH + M 

6.00E16 

-0.6 

0. 

8 

H + H = H 2 + M 

6.40E17 

-1. 

0. 

9 

H + 0 2 = H0 2 + M 

1.70E15 

0. 

-1000. 

10 

H0 2 +H = h 2 + 0 2 

1.30E13 

0. 

0. 

11 

H0 2 + H = OH + OH 

1.40E14 

0. 

1080. 

12 

H0 2 + 0 = OH + 0 2 

1.50E13 

0. 

950. 

13 

H0 2 + OH = H 2 0 + 0 2 

8.00E12 

0. 

0. 

14 

H0 2 + H0 2 = H 2 0 2 + 0 2 

2.00E12 

0. 

0. 

15 

H + H 2 O 2 = H 2 *f HO 2 

1.40E12 

0. 

3600. 

16 

0 + H 2 0 2 = OH +H0 2 

1.40E13 

0. 

6400. 

17 

OH + H 2 0 2 = H 2 0 + H0 2 

6.10E12 

0. 

1430. 

18 

M + H 2 0 2 = 20H + M 

1.20E17 

0. 

45500. 

19 

0 +0 = 0 2 +M 

6.00E13 

0. 

-1000. 

20 

N + N = N 2 + M 

2.80E17 

-0.75 

0. 

21 

N + 0 2 = NO + 0 

6.40E9 

1.0 

6300. 

22 

N + NO = N 2 + 0 

1.60E13 

0. 

0. 

23 

N + OH = NO + H 

6.30E11 

0.5 

0. 

24 

H + NO = HNO + M 

5.40E15 

0. 

-600. 

25 

H + HNO = NO + H 2 

4.80E12 

0. 

0. 

26 

O + HNO = NO + OH 

5.00E11 

0.5 

0. 

27 

OH + HNO = NO + H 2 0 

3.60E13 

0. 

0. 

28 

H0 2 + HNO = NO + H 2 0 2 

2.00E12 

0. 

0. 

29 

H0 2 + NO = N0 2 + OH 

3.43E12 

0. 

-260. 

30 

H + N0 2 = NO + OH 

3.50E14 

0. 

1500. 

31 

0 + N0 2 = NO + 0 2 

1.00E13 

0. 

600. 

32 

M + N0 2 = NO +0 

1.16E16 

0 . 

66000. 
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Table 4.2 Third Body Efficiencies 


Reaction 

Third Body Efficiencies 1 

H + OH + M = H 2 0 + M 

h 2 

1.0 

h 2 o 

6.0 

H + 0 M = OH + M 

h 2 

1.0 

h 2 o 

5.0 

h + h + m = h 2 + m 

h 2 

1.0 

h 2 o 

6.0 

H + 0 2 + M = H0 2 + M 

h 2 

2.0 

h 2 o 

16.0 

M + H 2 0 2 = 20H + M 

h 2 

1.0 

h 2 o 

15.0 


Table 4.3 9-Species and 18-Reactions Hydrogen-Air Model 


Species are H 2 , 0 2 , H 2 0, OH, H, O, H0 2 , H 2 0 2 and N 2 . 


j 

Reaction 

1 

H 2 + 0 2 = OH + OH 

2 

OH + H 2 =H 2 0 + H 

3 

H + 0 2 = OH +0 

4 

O + H 2 = OH + H 

5 

OH + OH = H 2 0 + O 

6 

H + OH = H 2 0 + M 

7 

H + H= H 2 + M 

8 

H + 0 2 = H0 2 + M 

9 

H0 2 +H = h 2 + 0 2 

10 

H0 2 + H = OH + OH 

11 

ho 2 + O = OH + 0 2 

12 

H0 2 + H 2 = H 2 0 2 + H 

13 

ho 2 + OH = H 2 0 + 0 2 

14 

ho 2 + ho 2 = H 2 0 2 + 0 2 

15 

H + H 2 0 2 = H 2 0+ HO 

16 

O + H 2 0 2 = OH +H0 2 

17 

OH + H 2 0 2 = H 2 0 + H0 2 

18 

M + H 2 0 2 = 20H + M 


All other third bodies have efficiency of 1 .0 
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Table 4.4 7-Species and 7-Reactions Hydrogen-Air Model 


The species are N 2 , O 2 , H 2 , OH, H 2 O, O, and H. 


j 

Reaction 

1 

0 2 + H 2 ^ OH + OH 

2 

0 2 + H ^ OH + O 

3 

H 2 + OH r=± H 2 0 + H 

4 

H 2 + 0 ^ OH +H 

5 

OH + OH ^ H 2 0 + 0 

6 

OH + H + M ** H 2 0 + M 

7 

h + h + m^h 2 + m 
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Chapter 5 

METHOD OF SOLUTION USING SHOCK CAPTURING 

This chapter deals with the solution procedure using shock capturing method. In this 
work we used both shock capturing method as well as shock fitting method. Shock-fitting 
will be discussed in the next chapter. In both the methods the standard MacCormick’s 
central difference algorithm was used. 

5.1 Shock Capturing Method 

In shock capturing the outer boundary of the coordinate system is outside the bow 
shock wave. Here, the shock comes naturally out of the finite-difference calculations, 
showing up as a rapid change of flow properties across several grid points. It is not 
treated explicitly as a discontinuity, and the oblique shock relations are not used. Once 
the thermodynamic properties, diffusion coefficients and chemical production rates have 
been defined, the governing equations can be solved numerically. The finite difference 
solution procedure is discussed in the next section. 

5.2 Finite-Difference Solution Method 

To solve the governing equations 3.1 with the finite-difference scheme, the equations 
must first be transformed from the physical domain (jc.y) in which they are written to 
an appropriate computational domain (£,??). The grid is kept uniform in both x and y 
directions. 

To transform the governing equations from the physical to the computational domain, 
fluxes F and G are first written in functional form and differentiated with respect to the 
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computational coordinates. Therefore, given F = F(x,y) and G = G(x,y), and we have 


Ff - F x x t +F y y ( 


Fy — Fjr^+FyJ/^ 

Then, substituting F y from Eq. (5.2) into Eq. (5.1) and simplifying gives 

r. F^-FvVt 

* X T 


where 


J = x (Vv - Vi x v 


(5.1) 

(5.2) 


(5.3) 


(5.4) 


is the Jacobian of the transformation. 

Proceeding in the same manner for G gives 

n G V X ( —G^X V 

G„- - 


(5.5) 


Finally, substituting Eqs. (5.3) and (5.5) into Eq. (3.1) gives the governing equations in 
the computational domain 


dU dF dG ~ 
dt dr} 


(5.6) 


where, 

U = JU, H = JH 


F = y„F - x„G (5.7) 

G = x ( G - y ( F (5.8) 


Here (x* x n> y ( , y n ) are the transformations metrics that form the inverse Jacobian metrics, 
and J is the Jacobian of the transformation. The metrics can be computed numerically 
once the physical coordinate grid has been prescribed. 


53 



5.3 MacCormack’s Implicit (Partial) Method 

The MacCormack’s [48] finite-difference method is a predictor-corrector scheme of 
the Lax Wendroff type. The scheme is second-order accurate in time and space, which 
results in a spatially and temporally discrete simultaneous system of equations at each grid 
points. To solve the governing chemically reacting flows equations, the spatial derivatives 
must first be discretized, and then an approximate temporal discretization must be chosen 
to advance the equation ahead in time. The temporal scheme must be chosen carefully 
because the system of partial differential equations describing chemically reacting flows 
can be stiff because of the highly disparate time scales that exist among the equations. 

The governing equations can be stiff because of the kinetic source terms contained 
in the vector H and because of the diffusive terms in the vectors F and G. Only the 
kinetic terms introduce stiffness in this work; spatial stiffness is controlled by the choice 
of the grid. To deal with the stiff system, the kinetic term is computed implicitly. Since 
it is only partial implicit so only the U vector and source term vector H are advanced in 
time level n+ 1. Thus in temporally discrete form we have 


U n+1 - U n • 

_±j 

At 



(5.9) 


Therefore, 


u r,r = u "> - ai 





(5.10) 


Equation (5. 10) above is non-linear and so a linearization procedure must be followed. 
A linearization procedure based on Taylor series expansion is used. Thus H is linearized 
by expanding it in a Taylor series in time. All flux vectors are expressed in terms of the 
flux vector U. Thus H is linearized by Taylor series as 


H;,r= H” J+ 


m 

dt 


At + i?(At) 2 


(5.11) 
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In order to rewrite 32- in terms of the gradients of the flux vector U we know that 

H = /(U,6,£*,U 


The chain rule of differentiation yields the following relation 

dH _ dHdU dHdh dHo dH djy 

dt QV dt d(t dt d( x dt d£ y dt 


(5.12) 


For many applications the grid is independent of time. Therefore the above equation 
reduces to 


m _ m du 

dt ~ <9u dt 

Therefore, substituting this in Eq. (5.11) gives 

>*s. 

H?} 1 = H" j + + rf(A f) 2 

, ’ J dU dt v ’ 


(5.13) 


(5.14) 


The partial derivative is approximated by a first order forward difference expression as 


d\J U n+1 — U" n/A . AU n/A , 
Ht = Ai + = at + <(Ai) 


(5.15) 


Substituting the preceding Eq. (5.15) in Eq. (5.14) gives 


H S3 1 = HJ,- + — AU + d(A ty 


(5.16) 


Terms like are called Jacobian metrics. In Eq. (5.15) above, the term AU n+1 = 
U n+1 — U n is the change in flow properties per time step. The finite difference equations 
shall be formulated in terms of AU which is referred to as the delta formulation. 


Next the linearized Eq. (5.16) is substituted in Eq. (5.10) to yield 


(5.17) 


This gives 


AU dF\ dQ\ ~ 3H att 

aT+ 1*) + (w) =H - + au AU 

\ / i ,j \ / i,j 


(5.18) 



Let (||j) n = K n equals the Jacobian metric of H with respect to U. 

Als0 (^), J + — H "j = R" ; equals the steady state residuals. Therefore Eq. 

(5.18) reduces to 

AU n+1 

=-A — + R" - K n AU n+1 = 0 

At 

AU” +1 + AtR n - K n AtAU n+1 = 0 

AU" +1 [I - At K n \ = -AtR n (5. 19) 

where I is the identity metrix. Once the temporal discretization has been used to 
construct Eq. (5.19), the resulting equations are spatially differenced by using unsplit 
MacCormack’s scheme. The modified MacCormack’s technique becomes 

53.1. Predictor Step 

The predictor step is given by 

[I - AtK*j] AU^f = -AtR”j (5.20) 

where, R”j represents a forward spatial difference of R. Also AU”j 1 = U"* 1 — 
where n + 1 is the temporary advanced level of time. For the second step (corrector) a 
backward differencing for spatial is used for R. 

53.2. Corrector Step 

The corrector step is given by 

I - AtKff'] AU" j 1 = -AtR^f (5.21) 

where this time a backward spatial differencing for R is used. Because, 

U n+1 - U" = AU n+1 


Therefore, 


AU n+1 - -AU n+1 = -AU n+1 
2 2 
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This gives. 


or. 


-AU" +1 = -AU n+1 
2 2 


AU n+1 = AU” +I 



since u), remains nearly constant over a time step. For accuracy, it is required that the 
chemical time step be chosen such that no change in mass fraction greater than 0.01 
occurs over that time step. The computational time step At is then chosen to be the 
minimum over all grid points of the fluid and chemical time steps. 
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5.4 Boundary and Initial Conditions 

The governing Eqs. (3.1) require boundary conditions along all the three boundaries. 
For the problem to be considered the inflow boundary is always supersonic, so the 
velocity, static temperature and pressure, and species are specified and fixed there. At 
the supersonic outflow boundary, all flow quantities are extrapolated from interior grid 
points. Although full Navier-Stokes (N-S) equations are used, the slip conditions are used 
to numerically simulate the inviscid flow. A flow tangency or slip boundary condition 
is implied on solid wall. The wall temperature and pressure are extrapolated from 
interior grid points. Initial conditions are obtained by specifying free-stream conditions 
throughout the flow field. 

5.5 Artificial Viscosity 

The Lax-Wendroff type schemes are inherently unstable and, hence, higher-order 
numerical dissipation terms are often necessary to get a stable solution. For a non-reacting 
flow field, an artificial viscosity based on temperature and /or pressure is traditionally 
used, but in chemically reacting flows, in addition to temperature and pressure gradients, 
one can also have very strong species concentration gradients. To suppress the numerical 
oscillations in the induction zone where the gradients in the concentration of reactants 
and products are very strong, additional artificial viscosity based on water mass fraction 
is used similar to the one used by Singh et al. [35] 
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Chapter 6 

METHOD OF SOLUTION USING SHOCK-FITTING TECHNIQUE 

The bow shock generated by an object in a supersonic/hypersonic flow field may be 
selected as the outer boundary of the domain and determined as a part of the overall 
solution. This procedure is known as the shock-fitting method. In problems like 
shock-induced combustion where physical instabilities are present, the shock capturing 
methods will smear some of the instabilities. Thus shock-capturing methods, when used 
in complicated problems of practical interest, would not reproduce many intricate flow 
features. On the other hand, in the shock-fitting approach, one knows the precise location 
of the discontinuity which acts as one of the boundaries of the flow field, across which 
Rankine-Hugoniot conditions are applied. This approach avoids taking differences across 
the shock and the smearing of the shock that occurs in the shock-capturing method. 
There are some obvious advantages of shock fitting over shock capturing. Shock fitting 
requires far less grid points compared to shock capturing. In shock capturing the bow 
shock becomes a smeared shock surface and, requires more grid points for the extension 
of the grid in the free-stream region. This adds to the savings in computational time 
in shock fitting as compared to shock capturing. Since very small dissipation is needed 
in shock fitting, the intricate details of the flow can be reproduced, as the dissipation 
does not smear the important flow features. This chapter deals in details the solution 
procedure used with shock- fitting method. 

6.1 Time-Dependent Navier-Stokes Equations in Strong Conservation Law Form 

Navier Stokes equations in vector form given by Eq. (3.1) and repeated here for 


59 



convenience is 


U t + F x + G V = H (6.1) 

The governing equations of motion are transformed from physical space (jc,y) to the 
computational space (£,17) by the following relations: 


T = t 

(6.2) 

\KS> 

II 

(6.3) 

'kT 

II 

(6.4) 


Applying the chain rule of partial differential equation gives the following expressions 
for the cartessian derivatives. 


d d d d 

dt ~ dr +m dr ) 

d_ _ d_ d_ 
dx - 

d_ _ d_ d_ 
dy ~ ^ y dt +T}y dr) 

Therefore Eq. (6.1) reduces to 


(6.5) 

( 6 . 6 ) 
(6.7) 


u r + + ritU v + £ r F{ + %F v + £yG{ + r) y G v = H (6.8) 

6.2 Computation of Metrics £t,Vt,Zx,Vx,£y,Vv 
From Eqs. (6.5)-(6.7) it is obvious that the value of the metrics (t,Vt^x^Vx^(yiVy 
must be provided in some fashion. In most cases the analytical determination of the 
metrics is not possible and must be computed numerically. Since the step sizes in 
the computational domain are equally spaced x^, x v etc can be computed by various 
finite difference approximation. Thus, if the metrics appearing in Eqs. (6.3)-(6.5) can 
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be expressed in terms of these derivatives, the numerical computation of metrics is 
completed. To obtain such relation, the following differential expression is considered 


, dt dt.^dt 

6 ‘ ~ a7 + 3{ d{ W ? 


But according to Eq. (6.2) 


dt_ 

dr 


= 1 


and 


dt _ dt _ 

d~t~d^~ 


Therefore, 


dt = dr 


dx = x r dr + x^d£ + x v drj 
d y = y T dr + y f d£ + y v dr) 
Writing Eqs. (6.9)-(6.11) in matrix form we have 


dt 


‘ 1 0 O' 


dr 

dx 

= 

X 7" X £ X fj 


dt, 

dy _ 


J lr Vt Vf) . 


dr] 


Reversing the role of independent variables, we can write 

dr = dt 

d£ = £ t df + £ z dx + £ y dy 
dr] = T] t dt + r] x dx + r] y dy 

Therefore, 


dr 


i 

o 

0 

r-H 

1 


dt 

dt 

= 

tt tx ty 


dx 

dr] 


_Vt Vx Vy. 


dy . 


(6.9) 

( 6 . 10 ) 

( 6 . 11 ) 

( 6 . 12 ) 


(6.13) 

(6.14) 

(6.15) 

(6.16) 
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Comparing Eqs. (6.12) and (6.16) we have 


"1 

0 

0 ' 


' 1 

0 

o ‘ 

6 

£x 

£y 

— 

X T 

H 

x v 

jit 

Vx 

Vy. 


JJr 

Vi 

Vv. 


w 

det | denominator | 


where 


Let 


Therefore, 


[X] = 


x i 

x v 


0 

0 


0 

0 ' 

- 


Vn, 


y* 




Xjf 


X T 

x v 


i 

0 


l 

0 ' 


Vr 

L 

y*. 


Vr 

Vv 


x T 

l 

x v. 


X T 

x t 


i 

0 


1 

0 


Vr 

V(. 


.Vr 

yt. 



x t. 

- 


J = 


1 

det |denominator| 


£x — Jj/»J 

£y = 


£ ,t — J {^tUt) jVt) 


Vx = -Jy* 


Vy = 


T)t = j {x r ys - X{y T ) 


In the evaluation of Eq. (6.8) many terms which adds to zero analytically but numerically 
are not zero (known as GCL terms) have been neglected. Now we evaluate those terms 
and add it to the source term H. We can write 


Ut 


dU , dU d\J 
8 t + ** d( + m dr, 


= + (tue.it - u(e.)t} + {[u„]„ - u<„),} 
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Also, 


I\r — + T] jF, 


= (F6) { - F({,) 4 + (F,,)„ - Ffo,), 
= (F&)« + (Fi).), - [F(&) { + F(,,), 


similarly. 


iy = (G(,) e + G(%), - G[((,) { + 


Substituting the values of Uj, F x , and G y in Eq. (6.1) gives 

+ [U6 + F£ x + + [U Tj t + Ft) x + GT] y ] r 

GCL TERMS 


(6.17) 


= H + [U(f,) { + U(„),] + [F({,) { + Ffo,),] + [G(£,) { + G(,,)J 
Dividing Eq. (6.17) through the Jacobian of the transformation and substituting the values 
of £x,( y ,(t, Tfr, and rjt and rearranging terms we have 


(JL\ 

( * 

\ 


( G 

\ 


[{-x T y t, + Xr,y T } U + y v F 

— x y G] I 

+ 


- y^F + G] 

V Jr 


J 



/ 



H 





GCL TERMS 

^ - 


= J lR 4- [u{-z r y„ + x v y r } i + U{x r y f - 2 *y r } J + [F(y„) f - F(y^ - G(a^ + G(x^) f? 

(6.18) 


Therefore, the final form of the governing equation in the computational domain with 
time-dependent terms is given by 


dU &F dG > 
dt + d£ + d v ~ H 


In Eq. (6.18), the terms that add to zero analytically but numerically are not zero are 
referred to as geometric conservation law (GCL) correction terms. 
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6.3 Shock Conditions 


The flow conditions along the leftmost boundary are those conditions which exist 
immediately downstream of the bow shock as determined by Rankine-Hugoniot relations. 
Therefore, this boundary is allowed to move with the bow shock as the later moves 
towards the steady state position. Figure 6.1 shows the coordinate transformation used 
in shock fitting procedure in which the following notations are employed. 

Let 6 equals the angle the radius makes with the horizontal (i.e. body angle). 
a equals the angle the radius makes with respect to body (i.e. tangent to the shock 
makes w.r.t. tangent to the body). 

Taking projections of unit vectors i r , io on n and s we have. 


s = ig cosa + i r sina 
n = i r cosa — ig sina 

Let A0 is the small increment in the body angle. 


tana = 


cosa = 


sina = 


_ Aii - JL 


r s A0 


= — X 


do 


l n 

r 9 





\ A +( t 5-) 2 

Let V\ represent the vector component of the fluid velocity normal to and measured 
with respect to the moving shock. 

Therefore, one may express the following: 


Now i r ■ n = 


Voo -Us 


n 


}n = {(««,*; -n + ««,*■; -n-tf.-n) • n} (6.19) 


n cosa 


Let 



( 6 . 20 ) 
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J2. , 

Therefore, sina = and cosa = ^ 


l T ■ Tl = 


0 


and 

— * —4 

iff ■ n = - iff 

Equation (6.19) reduces to 

V\n = 


n\ cos (90 — a) = — sin a = 


r _seh 

0 


v oo r se/ r s TT fj ~ 

T~~r 


Tl 


( 6 . 21 ) 


( 6 . 22 ) 


Let r St equals the radial shock velocity. 

Therefore, magnitude of the shock velocity in the direction normal to the body (i.e. 
along the radius) is given by 

dr s U s 


r St = IT = — = U *P 
ot cos o 


Hence, Eq. (6.22) reduces to 

r„ = Vi 


l + (-^ 

r 3 


n i/2 


+ Voq 


"”(r?) 


(6.23) 


The derivative that appears in Eq. (6.23) is evaluated by using the second-order 
central difference formula 

. Ktn-'VJ 


■8 0 / 


(2 < j < nny - 1) 


'0) 2A 9 

At the beginning of the predictor step, the shock wave radial distance is computed from 
the Euler predictor equation 

dr s 

r S( t +At) = r S(t) + 


or. 


r s +1 = r n s + At r« 


(6.24) 
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Therefore, 


= VJ 


Vi = Vi ■ n 


i r cosa — ig sin a 


= [Vicosa]i r + [— Visinajtfl 


V^Zjr "j” U}Z 


where, v\ equals the component of fluid velocity Vi normal to the body ( i.e . along i T 
direction), and ui equals the component of the fluid velocity V\ tangential to the body 
(i.e. along iff) direction. Therefore, 


ui = — 


fr) 




T s t v oo + ttoo 


V r s ) 


and 


Vj = Vi cos a — 


St 


-»„ +Uoo (^)]-r= 


, i+ (^) 

Let u is equal the velocity component tangent to the shock (i.e in s direction) and v in 
equals the velocity component normal to the shock (i.e. in n direction). 

Then by applying the shock jump conditions we have 


Pi^in — p2^2Tl 


Because tangential velocity is preserved, therefore 


U 2 s — 

Let V 2 equals the resultant velocity after the shock with respect to shock coordinates. 
Therefore, V 2 = v 2n n + u 2S s 

The component of V 2 along ig (i.e tangent to the body) is given by 

u 2 = V 2 ■ ig = v 2n n ■ ig + u 2S s ■ ig = — u 2n sino; + u 2S cosa (6.25) 
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Similarly, the component of V 2 along i r (i.e normal to the body) is given by 


V 2 = V2 • *r = (v 2n n + u 2S s) • ir = v 3Tt cos a + ursine* 


(6.26) 


^15 — t>ooir “t“ ^00^0 ^5 * ^ 


= Voo i r ■ s + iiooiff ■ s — u s ■ s 


— Voo sin a + u 0 0 cos a 


(Because U s • s = 0; i T ■ s = sin a; ig ■ s = cos a ) 


Therefore, 


U ls — U2S — 


U °° ("rf ) + U ° 


'1+ Tf- 


From Eq. (6.25) we have. 


(6.27) 


U2 = i> in sina + u 13 cosa 

P2 

Pi . \Voc 7 frr r *. , COSQ L f ra A_LTr 

P 2 0 0 0 \ 0 l \r,/ J 

Rearranging terms we have, 

Tangential component of the fluid velocity after the shock with respect to body is given 


( Pl \ r *. - 
U 2 = U 00 -[l-^-jx i jr- 


- (S)1 


(6.28) 


Similarly, normal component of fluid velocity after the shock with respect to body co- 
ordinates is given by 


v 2 = v 2n cos a + u 2S sin a 
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Substituting values and rearranging terms, we have 



(6.29) 


Pressure behind the shock is obtained from the MacCormick’s scheme at the predictor 
level. Once the pressure is known behind the shock, the normal component of the flow 
velocity ahead of the shock (measured w.r.t. shock) can be related to the pressure behind 
the shock by manipulating the oblique shock relations, which are 


P i«i = P 2 V 2 


Pi + Pl^j 2 = P 2 + p 2 A 


h 1 


u, = Uo 


V 1 l V 2 

+ T ~ 1,2 + T 


where Vi and V 2 are resultant velocities. 
From Eq. (6.31) we have 


p.vf - p?vl = p 2 - p, 


(6.30) 

(6.31) 

(6.32) 

(6.33) 


From Eq. (6.30) we have 


Pl_ 

P2 V, 


P2 - Pi 



fhvl 



= (/>* 



(6.34) 


Therefore, 


2 _ ( P2 ~ Pi \ / P2\ 
* \P2 Pi / \Pl/ 


Pi 

Pi 


21 — i 
Pi 


£2 

pi 



- - 1 (— \ - + 1 — + 
2 \Pl/L7-lpi 
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Multiplying numerator and denominator by we have. 




P2 7 - 1 


2 /PiLpi 7 + !j 


Similarly, 


Therefore, 


2 = £ 1z pi/M 

P2~Pl\P2 ) 

V 2 V 2 

*' + t=‘ 2+ f 


c,(r 2 -Ti) = i(v , 2 - v?) 


Velocity triangles upstream and downstream of the shock gives 

V\ ~ v 2 = W + - ( v 2 + u a) = v i - 

Equation of state gives. 


'Y a 

ft = c p T = -t-rRT = 


7 P 


7 — 1 7 — 1 7 — 1 p 

Substituting Eqs. (6.37) and (6.38) into Eq. (6.35) we have. 


7 


( El 

_ pA 

- ^( u i ~ v \) = \ 

>2 -pii 

El 

pi' 

\p 2 

pi) 

2 2 

to 

1 

.Pi 

P2. 


1 (P2 ~Pl)(p2 

2 PlP2 


After rearranging we have, 

27[plP2 - P2P1] = (7 - 1)(P2 -Pl)(P2 - Pi) 

P 2{(7 ~ 1 )P 2 - (7 + l)Pl} = Pi {(7 - l)Pl - (7 + 1 )P 2 } 


P 2 
Pi 


I±l£2. _ 1 

7-1 Pi 1 

2±1 _ £1 
7-1 Pi 


Equation (6.39) can also be written as 

p 2 { — 27P1 - (7- 1)(P2 — Pi)} = Pi {(7 - !)(p 2 -pi)- 27 p 2 } 


(6.35) 


(6.36) 


(6.37) 


(6.38) 


Pi) 


(6.39) 
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(6.40) 


P2 

Pi 

P 2 _ 
P 1 


(7 + 1 )P 2 + (7 ~ l)Pl 
(7 + l)Pl +( 7 ~ 1 )P 2 


2±lEl + I 
7-1 Pi Z 

2±I _l El 
7-1 Pi 


El _i_ IzI 
Pi _ 7+1 

1 4. l^i£2. 
T 7+1 Pi 


Equations (6.24), (6.28), (6.29), (6.35) and (6.40) when written in the notations of the 


advanced time level in terms of the preceding time level can be written as 


r ^+l 

1 3 


= r; + At r» 


(6.41) 


yn+l _ 
nnx,j 


{ 7 ~l~ 1 A Poo PnnxJ + 7 1 

\ 2 / Poo Poo 7 " 4 * 1 


~n+l 

J nnxj 


EjULLlL 1 

Poo ~ 



poo 


1 + 


Ini 

7+1 Poo 


= Uqo ~ 

nnx,j 00 


Poo 

n-f-T 

nni,; 



?; n + l _ 1 

u nm,; ” u oo 



Poo 

p n+1 

rnnx y j 



(6.42) 


(6.43) 


(6.44) 


(6.45) 


Note i is normal to the body and varies from i = 1 at the surface to / = wrt at the shock. 


Also j is along to the body and varies from j = 1 at the stagnation line to j = nny at 
the outflow boundary. 


6.4 Solution Procedure 


Solution procedures are followed in four steps as described below. 

STEP 1 : Initial Solution 

The initial conditions for this calculation are obtained by using an approximate 
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curve fit for the location and shape of the bow shock. Newtonian pressure distribution 
is used at the body. The approximate curve fit of Billig [46] is used to find r s and r se 
along the shock. To find the initial conditions immediately behind the shock, the radial 
shock velocity r St is set equal to zero, and Eqs. (6.19) and (6.42)-(6.45) are used. The 
initial flow conditions on the wall are obtained by using the known wall temperature in 
conjunction with the pressure from the Newtonian pressure distribution. The initial flow 
conditions at the interior grid points are obtained by assuming a linear variation between 
the flow conditions immediately behind the bow shock and the wall conditions. 

STEP 2 : Predictor Step 

At the beginning of the predictor step, the shock-wave radial distance is computed 
from Eq. (6.41). The pressure immediately behind the shock (pnnx,j) is computed with 
the MacCormack scheme 


U n+J . = U n 

^ nnxj nnx j 


( pin l?" \ 

^ \ x nnxj+l XnnxJ) 


At_ 

Ax 




nnx ,j 


(6.46) 


After the pressure behind the bow shock is obtained, V"+ x j 7 and can be computed 

from normal shock relations given by Eqs. (6.42) and (6.43). Similarly, the components 
of the fluid velocity behind the bow shock can be found from the oblique-shock relations 
given by Eqs. (6.44) and (6.45). The remaining unknown T”^- is calculated with the 
equation of state. This completes the predictor step. 

STEP 3 : Corrector Step 

The corrector step is similar to the predictor step, except that the shock-wave radial 
distance is evaluated with the modified Euler corrector, which is 


r n + 1 — r « 
s O) S U) 


+ 


— (r n 

2 \ St U) 


+ r 


n+1 

'Hi), 


(6.47) 
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and Eq. (6.46) is replaced by the MacCormack corrector scheme, in which the usual 


backward difference for dG/dy is replaced by a forward difference given by 

u n +^ — 1 tj” . _L u n +^ — ( 

nnx,j 2 nnxj' nn x,j Ay\ nnx 'J r nnx,j-l) I 

- \ [£ ( G S - G ^_, ,) + 


(6.48) 


This completes the corrector step. 

STEP 4 : Boundary Conditions 

When the “shock-fitting” method is employed, the flow conditions at the outer compu- 
tational boundary are those conditions which exist immediately downstream of the bow 
shock as determined by the Rankine-Hugoniot relations. Consequently, it is necessary to 
include logic which will permit this boundary to move with the bow shock as the latter 
moves toward its “steady-state” position. After the calculation of the boundary condition 
at i = nnx (i.e., after the shock) has been performed by the shock-fitting method, the 
predictor or the corrector steps are initiated at the interior grid points. All other boundary 
conditions are calculated after the predictor or corrector step has been completed at all 
interior grid points. 

The flow conditions along the supersonic outflow boundary (i.e., at j = nny) are 
determined by using a second-order extrapolation of the interior grid-point data. 

Along the body surface, no-slip, zero-pressure-gradient, adiabatic, and noncatalytic 
wall boundary conditions were used. 

6.5 Artificial Viscosity 

It is well known that central finite-difference schemes for non-linear convection 
problems require the addition of artificial dissipation for stability and uniqueness. For 
problems with boundaries, the boundary conditions for the dissipative operators are not 
always obvious. In general case when no information about the field outside the boundary 
is available the dissipative operators used must be supplied with suitable boundary 
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Figure 6.2 Notation used for artificial viscosity. 


conditions. Such boundary conditions can be implemented either by extrapolating the field 
from interior grid points to the exterior grid points needed by the dissipative operators or 
by modifying the operators at the boundary so that no exterior grid points are addressed. 
For each extrapolation rule there is a corresponding boundary-modified operator and it is 
merely a matter of taste which representation is preferred. What is important is that the 
overall dissipative operator should be well behaved in the sense that it provides as much 
damping as possible without introducing unacceptable errors. 

In most applications where artificial dissipation is used, the dissipative operators are 
based on either 2nd-differences or 4th-differences in the various coordinate directions. 
They are of either constant coefficients or variable coefficients type. Since in shock-fitting 
we need very little artificial dissipation, generally fourth order damping is sufficient. 
Fourth-order explicit dissipation term of the following form is used in this work 

= [(U? +w - 4U” +U + 6U?, - 4U?_ U + U?_ 2J )] 

[(U&+2 - 4U" J+1 + 6U”, - + U?j_ 3 )] 

The negative sign is required in front of the fourth-derivatives in order to produce positive 
damping. The smoothening coefficient e e should be less than ^ for stability. The fourth 
derivative terms are evaluated using the finite-difference approximations shown on the 
right hand side of the above equation. The dissipative operators used are of the type 

DV itj = DiXJij + DjVij 

where A and D } are the corresponding 1-D operators in i and j respectively. Referring to 
Fig. 6.2 for various notations we can write the following fourth order difference equation 
in the interior domain for the radial direction (i direction). 

AU,, ,■ = -ee(Ui_ 2 j - 4U t _ w + 6U - 4U t+lj + U i+2 ,,) 
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This we do from i = 3 to i = nnx—2. The points remained to be solved in the i direction 
are at i = 1, 2 and i = nnx—l, nnx. Hence, 

AU lii = -e.[U 3li -2U 2l ,- + U lii ] 

AU 2j = -e e [U 4j - 4U 3iJ + 5U 2)J - 2Uij] 

•C^iUnni— 1 j = Ce[ — 2U nrU ;j + 5U nnl _ij — 4U nnl _ 2) j; + Unn*— 3,j] 

f^iUnnx,j = c e[U n nx,j — 2U nnl _i ) j + U nni _ 2,j] 

Again in the j direction we have the following fourth-order differencing used in the 
interior domain. 


DjUij — e e(U>,j-2 4Ujj_i + 6U,J — 4U, J+ 1 + Ui )j4 . 2 ) 

This we do from j = 3 to j = tiny— 2 and i = 1, nnx It needs to find out U t J at j = 1, 
2 and j = nny— 1, nny Therefore, 


= — e e [U,3 — 2U 1>2 + U t|1 ] 


Finally, 


= — €e[U,-,4 — 4U, i3 + 5U, i2 — 2U, ; i] 

^>U«,nny-l = ~ e e — 2U, )nn y + — 4Ui |T , n y_ 2 + Ui itm y _3 

Dj\} x ,nny — ~ e e[Uj jnn y — 2Ui >nn y_i + Ui >nn y_ 2 ] 

Ujj = U,j — e e [D,Ui,j + 


6.5 Validation of Code 

Before running the code with full chemistry and reactions, it was validated with two 
other existing sources. First it was validated with viscous shock layer code by Kumar and 
Graves [49] which also uses shock fitting method. Secondly, Bellig’s [46] correlations 
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for sphere were used to find the shock profile and the shock stand-off distance was 
compared with the results of the current code. The free stream conditions used for this 
validation were 


Moo = 5.11 

Poo = 42732 N/m 2 
Too = 292k 
7 = 1.4 

The computational mesh consisted of 101 equally spaced points in the radial i direction 
and 91 points in the tangential (J) direction. Figure 6.3 shows a plot of radial shock 
velocity (r st ) as a function of time step number n. These velocities should approach zero 
as the steady-state solution is approached. The solution moves rapidly toward a steady 
state initially with large time-dependent oscillations appearing. In the final part of the 
solution, the convergence is very slow and has a monotonic behavior. As the computations 
proceeds the residuals or L 0 c norm of the flow field based on density variation drops 
by 6 orders of magnitude and then stays constant. This is shown in Fig. 6.4. At this 
point the flow field has attained a pseudo-steady state, i.e. the flow field keeps oscillating 
about a mean. Figure 6.5 shows a plot of the pressure distribution around the wall of 
the blunt projectile non-dimensionalized with respect to the stagnation line pressure. The 
solid curve is the current numerical result. The circles are the fairing of VSL code due 
to Kumar and Graves [49]. Both the results are in excellent agreement. A comparison 
of the numerically predicted shock shape with the empirical result of Billig is shown in 
Fig. 6.6. The circular symbol represent the emperical results, while the numerical result 
is shown as a solid curve. Again the two results match perfectly well. Figure 6.7 shows 
the line plot for pressure along various j = constant lines and Figure 6.8 shows the line 
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plot for the density along various j = constant lines. The circles are the result with 
the present calculations and the triangles are the results with VSL code. The results are 
in good agreement with each other. Figures 6.9-6.11 show the pressure, density and 
temperature contours respectively with the present calculations on the upper half of the 
plot whereas the lower half is plotted with VSL code by Kumar and Graves [49]. Both 
the shock profile as well as the shock stand-off distance matches well with the VSL code. 
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Figure 6.3 Plot of radial shock velocity (r St ) vs time step number n. 
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Figure 6.4 Residual history of computations. 
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Figure 6.6 Comparison of numerically predicted shock shape with the emperical results of Billig. 
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Figure 6.9 Comparison of current numerical results for pressure (upper half) with VSL code (lower half). 
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Chapter 7 

RESULTS AND DISCUSSION USING SHOCK-CAPTURING METHOD 

In this chapter we shall be discussing the numerical simulations results for Lehr’s 
Mach 5.11 and 6.46 cases using shock-capturing method. Numerical simulations have 
been carried out for a 15 mm diameter spherical projectile and the following free-stream 
conditions: 

Moo =5.11 and 6.46 

Poo = 42663.2 N/m 2 (320 mm of Hg) 

Too = 292 K, M C j = 5.11 

The premixed fuel oxidizer mixture is taken as 2 H 2 + O 2 + 3.76N2. For the 
spherically blunt projectile, a fluid particle approaching the body first encounters the bow 
shock thus raising its temperature, pressure and density while lowering the velocity. This 
fluid particle then travels a distance equal to the induction length at elevated temperature 
before it encounters the reaction front. Once the ignition starts, chemical energy is 
released and another discontinuity known as reaction front is formed. In the induction 
zone, the temperature and the pressure remain relatively constant at the post shock 
conditions, while the concentrations of radicals build up very rapidly. In the stagnation 
zone, due to large residence time, it attains equilibrium while away from it, flow is in 
non-equilibrium. 

7.1 Mach 5.11 Case 

The assumption of flow-field symmetry about the stagnation line is invoked and, 
therefore, only one half of the flow field is calculated. Calculations have been carried 
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Figure 7.2 Contour plot of pressure for Mach 5.11 with shock-capturing, 




out for a grid with 197 points in the circumferential direction and 152 points in the 
normal direction. Figure 7.1 shows the typical grid which contains 197 x 152 grid points 
(197 normal to the body and 152 along the body). For clarity, every fourth grid point is 
shown in the figure. This grid was chosen based on the earlier work by Ahuja et al. [25], 
where the flow field was shown to be adequately resolved with this grid. For the present 
case of stoichiometric hydrogen-air mixture, the Chapman- Jouget velocity is the same as 
the velocity of the projectile for the Mach 5.11 case. Unsteady flow phenomenon can 
occur if the free-stream velocity of the projectile is around the C-J detonation velocity 
of the mixture. The residuals dropped by three orders in 12,000 iterations and then 
remained constant. 

Figure 7.2 shows the pressure contours. It shows the bow shock (similar to non- 
reacting flows) and the expansion waves as the flow expands over the body. The peak 
pressure occurs near the stagnation region where the shock is nearly normal and decreases 
away from it as the flow expands. The contours show some oscillations, which are caused 
by the fluctuating shock and reaction front and will be discussed later. The reaction front 
is not visible in this figure since the pressure stays constant across the reaction front. The 
shock standoff distance compares well with the Lehr’s shadowgraph (Fig. 2.2). 

The density contours (Fig. 7.3) show a very complex flow field with the presence of 
two discontinuities. The outer discontinuity is the bow shock and the inner discontinuity is 
the reaction front. These two fronts can be seen separated from each other by the induction 
distance. The separation (i.e., the induction distance) is minimum near the stagnation line 
and increased away from it. This is because near the stagnation line, bow shock is almost 
normal and, hence, the post shock temperature is maximum; thus, induction distance 
is minimum. Away from the stagnation line, the shock strength decreases, thereby 
decreasing the post-shock temperature and, hence, increasing the induction distance. 
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About one nose radius downstream, the induction length becomes infinite, i.e., infinite 
streamwise separation between shock and the reaction front. Therefore, it is not possible 
to form an oblique detonation wave under these circumstances. Also, the oscillations in 
the reaction front can be seen clearly. As the computations proceeds, the L <*, norm of 
the flow field based on density variation drops by two orders of magnitude and then stays 
constant. At this point the flow field has attained a psuedo-steady state, i.e., the flow 
field keeps oscillating about a mean. As noticed earlier, the density increases just after 
the shock (peak value) and then decreases as the flow passes through the reaction front. 

Figure 7.4 shows the temperature contours and the corresponding enlarged view. 
The bow shock and the reaction front can be seen clearly. Following a stream line into 
the stagnation zone, it is seen that temperature jumps across the shock and then stays 
constant in the induction zone. Past the induction zone, due to exothermic nature of H 2 - 
air reaction mechanisms, the temperature increases almost instantaneously reaching about 
1 1 times the free-stream value. The pulsation in the flow field can be vividly seen here. 

Figures 7.5 shows the water mass fraction contours. The reaction front can be seen 
in this figure. To better show the instabilities, an enlarged view is also shown. At the 
end of the combustion zone, the temperature is high enough to start the combustion. As 
the reaction proceeds, the water mass fraction increases rapidly. The oscillations similar 
to temperature and density profiles can be seen here. The instability is characterized 
by an almost regular periodic wave motion having a constant frequency and amplitude. 
Similar instability has been observed experimentally in Lehr’s work. The contour plots 
shown here show the spatial variation of instability at one point in time. It also shows 
that instabilities are not restricted to the reaction front only but are convected toward the 
body, thus affecting the entire flow field. Figure 7.6a shows the line plot for pressure 
along various j = constant grid lines. As the flow crosses the shock, it encounters the 
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pressure jump. The pressure decreases slightly after the shock. The Von Neumann spike, 
which is characteristic of reacting flows, is also visible. The post shock oscillations in 
pressure along stagnation line has also been observed in Ref. [14], Figure 7.6b shows 
the temperature along various j = constant grid lines. The post shock stagnation point 
temperature is 3150°K, which compares well with Ref. [15]. As the gas encounters 
the bow shock, the temperature increases abruptly. Immediately after the shock, the 
temperature stays constant for a short distance and then begins to increase due to 
exothermic reactions. The induction zone decreases with increasing temperature, as 
chemical energy release will be faster for higher temperatures. 

To help understand the temporal nature of these instabilities, attention is now 
focussed on the time history of physical variables along the stagnation line only. The 
presentation of results in this form also allows comparison of the numerical results with 
the wave-interaction model originally proposed by McVey et al. [8] and further modified 
by Matsuo et al. [18]. This is a one-dimensional model which is used to explain the 
instability mechanism. The model successfully explains the main features of the flow 
field as identified in the numerical simulations as well as in schlieren pictures. The 
essential features of this model are shown in Fig. 7.7, which also shows the x-t plot 
for the computed water mass fraction along with an overlay of pressure. Also, one 
can see the unsteadiness in the reaction front. This unsteadiness originates from the 
induction zone near the stagnation line and then travels downstream. First, the contact 
discontinuity approaches the original reaction front. The hot gases behind the contact 
discontinuity begin to react. The resulting pulse of the energy release front generates two 
compression or pressure waves, one of which propagates upstream towards the bow shock 
and another downstream towards the body. The compression wave which propagates 
upstream interacts with the bow shock and produces a contact discontinuity behind the 
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bow shock. The bow shock is stronger after the interaction and, therefore, the gas is 
hotter on the upstream side of the contact discontinuity. The hot gases behind the contact 
discontinuity reduce the induction time and create a new reaction front thus generating 
another set of compression waves. The contact discontinuity then reaches the position 
of the original reaction front, extinguishing the reaction at this point because no more 
unreacted gas exists there. This reduces the rate of energy release thereby generating 
rarefaction waves. The reaction front begins to recede because of the increased induction 
time of the colder fluid. The other compression wave, travelling towards the body, gets 
reflected from the body to the bow shock, and interacts with it at about the same time that 
the most recently created compression wave arrives at the bow shock. The compression 
wave and the reflected compression wave from the body interact with the bow shock, 
thus providing a possible mechanism for the creation of another contact discontinuity 
i.e. secondary contact discontinuity. The gases being hotter on the upstream side of the 
contact discontinuity, start burning, again generating compression waves and the cycle is 
repeated. Matsuo et al. [18] also emphasized the importance of considering the reflection 
of the compression wave from the body in their calculations. The compression waves 
reflecting from the blunt body may not necessarily be in phase with the compression 
waves created by the new energy release front. Thus, once these reflected waves interact, 
they cause the flow to be not exactly periodic; however, the pulsating energy release 
front could still be nearly periodic. In some instances the original compression wave 
and the reflected compression wave may not hit the bow shock at the same time, thereby 
distorting some of the structure of the bow shock as clearly seen in Fig. 7.7. 

If one observes these oscillations very closely, it is seen that in each cycle the water 
mass production rate, which is also a measure of energy release, at first increases and 
then eventually decreases to zero. This is the point of extinguishment of the reaction 
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front. The reaction almost comes to a standstill at this point. Since the new reaction 
front generated has high energy release (and, hence, high water mass production rate), it 
sends new sets of compression waves, which propagate both upstream and downstream, 
and the above cycle is repeated. 

To further investigate the unsteady nature of the flow field, a Fourier analysis of 
the flow field was conducted. For this, data at various sample stations along the j = 
61 grid line were stored for 30,000 iterations to get good temporal resolutions. The 
grid used was 101x78, and all calculations were time accurate. Figure 7.8 shows the 
amplitude versus frequency plot obtained by using a Fourier transform. The flow field 
spectrum is well resolved, and it clearly shows the fundamental frequency of 1.2 e+6 
Hz and a peak amplitude of 0.004. It also shows subharmonics and high-frequency 
numerical noise. Experimental fundamental frequency, as given in Ref. [11], is 1.96e+6 
Hz. The discrepancies between the experimental and the numerical value could be due to 
improper grid resolution. The calculations were repeated for a finer grid (131 x 101). The 
grid aspect ratio was kept the same in both the cases. Figure 7.9 shows the frequency 
spectrum for the flow field with finer grid. The sample stations have the same physical 
locations as in the previous case. The dominant frequency is 2.0e+6 Hz, and the amplitude 
is 0.004. This frequency is in close agreement with the experimental value of 1.96e+6 
Hz. The above calculations were repeated once again for another finer grid of 197 x 152. 
The grid aspect ratio was kept the same and the sample stations have the same physical 
locations as in the previous cases. Figure 7.10 shows the frequency spectrum for this grid. 
The dominant frequency now is 2.1e+6 Hz., and the amplitude is 0.004. Thus, refining 
the grid has not changed the frequency and therefore, the oscillations in the reaction front 
are physical. The frequency prediction is very sensitive to not only the grid but also the 
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Figure 7.3 Contour plot of density for Mach 5.11 with shock-capturing. 




COTOUR PLOT FOR TEMPERATURE 


ENLARGED VIEW 


0.015 


0.010 

y 


0.005 


0.000 



T K 



3154.70 

3006.85 

2858.99 

2711.14 

2563.28 

2415.43 

2267.58 

2119.72 

1971.87 

1824.01 

1676.16 

1528.30 

1380.45 

1232.60 

1084.74 

936.89 

789.03 

641.18 

493.32 

345.47 




0.005 0.010 

x 



Figure 7.4 Contour plot of temperature for Mach 5.11 with shock-capturing. 
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Figure 7.5 Contour plot of water mass fraction for Mach 5.11 and corresponding enlarged view. 
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Figure 7.7 x-t plot of water mass fraction and pressure along the stagnation line showing 

proposed instability process. 
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Figure 7.8 Temporal frequency spectrum of water mass fraction for Mach 5.11 for 101x78 grid. 
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Figure 7.9 Temporal frequency spectrum of water mass fraction for Mach 5.11 for 131x101 grid. 




Figure 7.10 Temporal frequency spectrum of water mass fraction for Mach 5.11 for 197x152 grid. 


chemical kinetics model. The frequency of oscillations depends very strongly on the 
ignition delay, which in turn depends on the rate constants. Wilson [15] in his investi- 
gation showed that, by changing the reaction rate constants for reaction (2) [Table 4.3] 
from Jachimowski [45] model to those recommended by Wamatz [47], the frequency of 
oscillations changed from 530 KHz to 820 KHz (55% change). This is because Jachi- 
mowski’s model gives shorter ignition delay at higher temperature than Wamatz’s model. 
Therefore, slight difference between the experimental and computed frequency can be at- 
tributed to the uncertainty in the reaction rate constants, among other factors such as grid 
resolution and the numerical damping. 

7.2 Mach 6.46 Case 

The results for the Mach 6.46 case will now be presented. For the present stoichio- 
metric hydrogen-air mixture, the C-J velocity is Mach 5.11. Thus for the Mach 6.46 
which is a superdetonative case, the projectile speed is significantly above the detona- 
tion velocity of the mixture. The numerical simulations is carried out for the following 
free- stream conditions: 

Moo = 6.46 
Poo = 42732N/m 2 
Too = 292K 

Figure 7.11 shows the contour plot of temperature for the Mach 6.46 case. The bow 
shock and the reaction front can be clearly seen in the figure. They are coupled with each 
other near the stagnation line and up to about 60 degrees from the nose at which point 
they start decoupling from each other by the induction distance. This occurs because 
bow shock is almost normal near the stagnation line and the post-shock temperature is 
maximum. For Mach 6.46, a very small induction distance occurs as a result of the post 
shock temperature remaining significantly high up to some distance near the stagnation 
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zone. Away from the stagnation line, the induction distance is increased as a result of 
decreasing shock strength and post-shock temperature. A comparison with Fig. 2.3 shows 
that all the flow features are very well captured. The temperature further increases as 
the reaction proceeds due to the exothermic nature of the reaction. The peak temperature 
occurs at the stagnation point. Figure 7.12 shows the contour plot for density at Mach 
6.46. The bow shock has a very crisp and smooth profile. The reaction front, which 
is smooth up to about 60-65 degree from the nose region, is wrinkled with very small 
amplitude waves downstream. The maximum density is seen to be just after the bow 
shock, and minimum density is after the reaction front. Figure 7.13 shows the pressure 
contours, the bow shock, and the expansion waves as the flow expands over the body. 
Peak pressure occurs near the stagnation region where the shock is nearly normal, and 
peak pressure decreases away from the stagnation region as the flow expands. Figure 
7.14a shows the line plot for pressure along various j = constant grid lines. The flow 
encounters the pressure jump as it crosses the shock. The pressure decreases slightly 
after the shock. The post-shock oscillation in pressure along the stagnation line is also 
observed. Figure 7.14b shows the temperature along various j = constant grid lines. The 
post-shock stagnation point temperature is 3550° K. The temperature increases abruptly, 
as the gas encounters the bow shock. Immediately after the shock, the temperature stays 
constant for a short distance and then begins to increase due to exothermic reactions. The 
induction zone decreases with increasing temperature, as chemical energy release will be 
faster for higher temperatures. The contour plots for water mass fraction is shown in 
Fig. 7.15a. The temperature is high enough to start the combustion, at the end of the 
induction zone. The water mass fraction increases rapidly as the reaction proceeds. An 
enlarged view of the oscillations is presented in Fig. 7.15b and it shows the macroscopic 
behavior. The reaction front shows a smooth profile near the stagnation region, but 
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a corrugated pattern with extremely small amplitudes is observed 60-70 degrees away 
from the nose region. 

A Fourier analysis of the flow field was conducted to further investigate the macro- 
scopic steady nature of the flow field. Data at various sample stations along the stagnation 
line were stored for 40,000 iterations to get good temporal resolutions. The grid used 
was 197 x 152, and all calculations were time accurate. Figure 7.16 shows the amplitude 
versus frequency plot obtained by using Fourier transform. The flow-field spectrum is 
well resolved, and it clearly shows the fundamental frequency of 2.67e+6 Hz and a peak 
amplitude of 0.001. Harmonics and very high-frequency numerical noise are also shown. 
Experimental fundamental frequency for Mach 6.46 is not available. 

7.3 Effect of Nose Radius on Stability of Reaction Front 

The key parameters for the onset of periodic unsteadiness have been identified as (1) 
induction time, (2) reaction rate constant, (3) activation energy, (4) heat release rate and 
(5) projectile nose radius. In this study we shall be discussing the effect of nose radius 
on the stability of the reaction front while keeping the first four parameters constant by 
choosing a particular reaction model and by fixing the free-stream Mach number. 

7.3.1 Mach 5.11 and Projectile Diameter of 2.5 mm 

The diameter of the projectile was reduced to 2.5 mm while keeping the same 
free-stream Mach number of 5.11. Other free-stream conditions were also kept the 
same. Figure 7.17 shows the density contours and corresponding enlarged view. It 
clearly shows a very smooth bow shock and reaction front. Thus, reducing the projectile 
diameter caused the instabilities to disappear. Figure 7.18a shows the x-t plot of water 
mass fraction along stagnation line for the Mach 5.11 for a 15 mm projectile diameter. 
Figure 7.18b shows the x-t plot of water mass fraction along stagnation line for Mach 
5.11, but with a projectile diameter of 2.5 mm. In the former case the reaction front 
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clearly shows periodic oscillations, whereas the later case shows a smooth reaction front. 
The instabilities for Mach numbers lower than the C-J Mach number are due to the 
ignition delay. What causes these instabilities to disappear for the same Mach number 
of 5.11 but lower projectile diameter can be explained now. Because the projectile 
speed is same in both the cases so the physical induction length are also the same. The 
difference is relative scale of the induction length and the shock stand-off distance. In the 
later case of small projectile diameter the physical length scale is small, reaction period 
is relatively large and temperature increase rather gradually within that small physical 
scale. Therefore, compression waves are created gradually rather abruptly. Thus the 
unsteadiness disappears for the small diameter projectile. 
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Figure 7.11 Contour plot of temperature for Mach 6.46 with corresponding enlarged view. 
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Figure 7.12 Contour plot of density for Mach 6.46 with corresponding enlarged view. 
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Figure 7.13 Contour plot of pressure for Mach 6.46. 
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Figure 7.15 Contour plots of water mass fraction for Mach 6.46 and corresponding enlarged view. 







10 6 10 7 10 8 


Frequency (Hz) 


Figure 7.16 Temporal frequency spectrum of water mass fraction for Mach 6.46 for 197x152 grid. 
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Figure 7.17 Contour plots of density for Mach 5.11 and projectile diameter 2.5 mm and corresponding enlarged view. 
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Figure 7.18 Time history plots of water mass fraction for (a) Projectile diameter 15 mm (Mach number 5.11) and 

(b) Projectile diameter 2.5 mm (Mach number 5.11). 





Chapter 8 

RESULTS AND DISCUSSION USING SHOCK-FITTING METHOD 

In this chapter we shall be discussing the numerical simulations of Lehr’s Mach 5.11 
and Mach 6.46 cases using shock- fitting technique. A comparison of shock capturing 
method and shock-fitting method will clearly show the advantages of shock-fitting over 
shock capturing. Later we shall simulate Ruegg and Dorsey’s experimental work using 
the same technique with emphasis to the large-disturbance regime. 

The numerical simulations were conducted to simulate Lehr’s [11, 38] experimental 
results. The physical and free-stream conditions used in the simulation were: 

Moo = 5.11 and 6.46 
Poo = 42732N/m 2 
Too = 292K 

8.1 Mach 5.11 Case 

For Mach 5.11 case, calculations were carried out on a grid of 101x101. Figure 
8.1 shows the pressure contour plots for the Mach 5.11 case and Fig. 8.2 shows the 
corresponding enlarged view. The complicated wave pattern seen in Fig. 8.2 can be 
viewed as made up of two types of compression waves. One type of compression wave 
originates from the reaction front while the other has been reflected from the projectile 
body. The reflected compression wave interacts with the original compression wave and, 
at the point of interaction, two new waves are generated. This reflection produces the 
observed cell structure. The compression wave which moves towards the bow shock, 
overtakes it and causes the bow shock to move forward. Thus the kinks appearing 
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on the bow shock are due to some of its structure being distorted by the compression 
waves. Figure 8.2 also reveals that these pulsations in reaction front are strong near 
the nose region and dissipate near the shoulder regions of the projectile. This fact is 
further supported by Fig. 8.3, which shows the numerical value of pressure contours 
along the complete body. From the pressure values given, it is shown that the pulsation 
of compression waves that originate near the stagnation line are the strongest, and as one 
moves towards the shoulder region of the projectile, the pressure is reduced to almost 
atmospheric pressure. 

Figure 8.4 shows the pressure distribution along the stagnation line. The pressure 
increases from free-stream pressure to 1.332 e+06 N/m 2 as the flow passes the normal 
shock along the stagnation line. The flow then encounters the pressure wave (see Fig. 
8.2) which further raises the pressure to 1.375 e+06 N/m 2 . The pressure then drops to 
1.342 e+06 N/m 2 as it passes through the expansion wave. This pattern repeats itself as 
the flow encounters a series of compression and expansion waves. 

To help understand the temporal nature of these instabilities, attention is now focused 
on the time history of physical variables along the stagnation line. Figure 8.5 shows the 
time history plot of water mass-fraction along the stagnation streamline. It shows two 
discontinuities. The outer discontinuity is the bow shock, which shows little kinks in the 
structure, and the inner discontinuity is the reaction front. The highly periodic oscillations 
in the reaction front that originate near the stagnation line and then spread downstream 
are clearly evident. The bow shock is at 0.009224 meter and the projectile surface is 
at 0.0075 meter (projectile surface is not shown in the figure) from the center of the 
blunt body. As seen from the figure, the frequency of oscillation (which is inverse of 
the time period) can be calculated directly from the plot. This frequency is 2.0 MHz, 
whereas the experimental frequency from Lehr’s ballistic data for Mach 5.11 is 1.98 
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MHz. Also, the amplitude of the oscillations of reaction front is 8.0 e-05 meters. Figure 
8.6 shows the water mass fraction contours. The reaction front can be seen in this figure. 
The instability is characterized by a regular periodic wave motion having a constant 
frequency and amplitude. The reaction is complete near the body where the maximum 
water mass production can be seen. 

Density contours are shown in Fig. 8.7. It clearly shows the presence of two 
discontinuities. The outer discontinuity is the bow shock and the inner discontinuity 
is the reaction front. These two fronts can be seen separated from each other. Moving 
downstream the induction length increases because of lowering of post-shock temperature. 
The oscillations in the reaction front can be seen clearly. It is also seen that the density 
increases just after the shock and then decreases as the flow passes through the reaction 
front, in agreement with the experimental shadowgraph. 

Figure 8.8 shows the temperature distribution along the stagnation streamline. 
Following a streamline into the stagnation zone, it is seen that the temperature jumps 
across the shock and then stays constant in the induction zone. Past the induction zone, 
due to the exothermic nature of the H 2 -air reaction mechanism, the temperature increases 
rapidly reaching almost 1 1 times (3150 K) the free-stream value. To further compare the 
experimental data with the numerical data, Fig. 8.9a shows the computed shadowgraph of 
the Mach 5.11 case for the complete projectile and Fig. 8.9b is the corresponding enlarged 
view. It is seen that the bow shock and the reaction front are separated from each other 
near the stagnation line, and this separation keeps increasing downstream of the stagnation 
line. This is what was observed experimentally also. Also, the bow shock is quite smooth 
with very little waviness but the reaction front clearly shows a periodic behavior. 

By means of time history plots, a comparison of the numerical results with the wave- 
interaction model originally proposed by McVey and Toong [8] and further modified by 
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Matsuo and Fujiwara [18], can be made. In order to understand how the wave interaction 
model fits with the numerical results, we shall have to consider Figs. (8.10-8.12) together. 

Figure 8.10 shows the time history plot for the pressure along the stagnation 
streamline. Figure 8.11 shows the time history plot for temperature along the stagnation 
streamline. By comparing the actual model shown in Fig. 8.12b with the x-t diagrams 
of pressure and density shown in Figs. 8.12a and 8.12b, it is demonstrated below that 
the model proposed by McVey and Toong fits very well with the present numerical 
calculations. 

As shown in Fig. 8.12b, a contact discontinuity first approaches the original 
reaction front. The gases are hot on the upstream side of the contact discontinuity 
and comparatively cold on the lower side, as clearly seen in Fig. 8.11. These hot gases 
behind the contact discontinuity begin to react, generating compression or pressure waves 
that propagate both upstream and downstream, as seen in Fig. 8.10. The compression 
wave which propagates upstream intersects with the bow shock and produces a contact 
discontinuity behind the bow shock (Figs. 8.12c and 8.11). The bow shock is stronger 
after the interaction and, therefore, the gas is hotter on the upstream side of the contact 
discontinuity. The hot gases behind the contact discontinuity reduce the induction time 
and create a new reaction front, thus generating another set of compression waves. At a 
somewhat later time, the contact discontinuity reaches the position of the original reaction 
front, extinguishing the reaction at this point because no more unreacted mixture exists 
there. The rate of energy release is effectively reduced, which generates rarefaction 
waves as shown in Fig. 8.10. The reaction front begins to recede because of increasing 
induction time of the colder fluid. The compression wave traveling towards the blunt 
body gets reflected from the body, travels back to the bow shock, and interacts with it at 
about the same time that the most recently created compression wave arrives at the bow 
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PRESSURE CONTOURS 



Figure 8.1 Pressure contours for Mach 5.11 over the complete blunt projectile. 
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Figure 8.2 Enlarged view of pressure contours for Mach 5.11. 




Figure 8.3 Pressure contours 
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Figure 8.4 Pressure distribution along stagnation streamline for Mach 5.11. 
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Figure 8.5 Time history plot of water mass 
fraction along stagnation streamline for Mach 5.11. 
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Figure 8.6 Enlarged view of water mass fraction contours for Mach 5.11. 
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Figure 8.7 Density contours for Mach 5.11. 




Figure 8.8 Temperature distribution along stagnation streamline for Mach 5.11. 
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Figure 8.9 Computed shadowgraphs for Mach 5.11 (a) Complete projectile and (b) Enlarged view. 
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Figure 8.10 Time history plot for pressure along stagnation streamline for Mach 5.11. 
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Figure 8.11 Time history plot for temperature along stagnation streamline for Mach 5.11 
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Figure 8.12 Side-by-side comparison of time history plots for Mach 5.11 and the wave interaction model: (a) Pressure 
contours from simulation, (b) McVey and Toong wave interaction model, and (c) Density contours from simulation. 






shock. The compression wave and the reflected compression wave from the body interact 
with the bow shock, thus providing a possible mechanism for the creation of another 
contact discontinuity, i.e„ secondary contact discontinuity. The gases, being hotter on 
the upstream side of the contact discontinuity, start burning again generating compression 
waves; the cycle is then repeated as shown in Figs. 8.10. 

8.2 Mach 6.46 Case 

The results for the Mach 6.46 case are now presented. Due to close coupling of bow- 
shock and reaction front (i.e., small induction distance) at high Mach numbers, a finer 
grid was needed to resolve the flow field. Therefore, for Mach 6.46, a grid of 201 x 151 
was used with 201 points in the circumferential direction and 151 points in the normal 
directions. This is a superdetonative case, i.e., the projectile velocity is higher than the 
C-J velocity of the mixture. Figure 8.13 shows the pressure contours as well as the wave 
pattern similar to Mach 5.1 1. When compared with Fig. 8.2, it is clear that the frequency 
of the compression waves moving towards the body and moving towards the bow shock 
is much higher. The bow shock and the reaction front are almost coupled with each other. 
Figure 8.14 shows the density contours. For this case, a very small induction distance 
occurs as a result of very high post-shock temperature. The features of the bow shock and 
reaction fronts appear to be the same like the simulations with shock-capturing. Figure 
8.15 shows the pressure distribution along the stagnation streamline. There is a jump in 
pressure after the bow shock and then the pressure drops when the pressure wave meets a 
rarefaction wave. It increases again when it encounters another compression wave. After 
the energy release front, there is another jump in pressure. This pressure wave oscillates 
between a high value (when it encounters a compression wave) to a low value (when it 
encounter a rarefaction wave). Also, when compared with Fig. 8.4 for Mach 5.11, we 
see that there are more numerous oscillations in pressure for the Mach 6.46 case because 
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PRESSURE CONTOURS 



Figure 8.13 Pressure contours for Mach 6.46. 
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Figure 8.14 Density contours for Mach 6.46. 




Figure 8.15 Pressure distribution along stagnation streamline for Mach 6.46. 





Figure 8.16 Temperature distribution along stagnation streamline for Mach 6.46. 
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Figure 8.17 Computed shadowgraphs for Mach 6.46: (a) Experimental scale, (b) Enlarged scale, and (c) Highly enlarged view. 
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Figure 8.18 Time history plot of water mass 
fraction along stagnation streamline for Mach 6.46. 
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Projectile Surface 



(a) Shock capturing (b) Shock-fitting 


Figure 8.19 Comparison of shock-capturing vs shock-fitting: Time history plots of pressure for Mach 5.11 
along stagnation streamline (a) with shock-capturing method and (b) with shock-fitting method. 








of higher frequency compression waves generated. Figure 8.16 shows the temperature 
distribution along the stagnation streamline. The stagnation point temperature is 3550 K. 
The temperature increases abruptly as the gas encounters the bow shock. Immediately 
after the shock, the temperature stays constant for a short distance and then begins to 
increase due to exothermic reaction. Figure 8.17 shows the computed shadowgraph for 
density. Figure 8.17b and 8.17c shows the corresponding enlarged view. The bow shock 
and the reaction front remain coupled with each other up to about 60-65 degrees from 
the stagnation streamline, as observed experimentally, and then start decoupling. Also, 
the reaction front shows slight oscillations of very low amplitude. Figure 8.18 shows the 
time history plot for water mass fraction. The bow shock is at 0.00884 meters from the 
center of the blunt body. The distance between the bow shock and the reaction front is 
very small. Also, as is clear from the figure, the amplitude of oscillations of the reaction 
front is 2.5 e-05 meters which is quite small as compared with the Mach 5.11 case. The 
frequency of oscillations can be computed directly from this figure and it is found to be 
2.85 MHz, which is comparable with the earlier work [26]. Thus, Mach 6.46 case is a 
very high frequency but very low amplitude phenomena. 

8.3 Comparison of Shock Capturing vs Shock-Fitting 

As mentioned earlier that shock- fitting requires very small dissipation which does 
not smear the important flow features and therefore the intricate details of the flow field 
can be reproduced. This is clearly shown in Fig. 8.19 which is the x-t plot for water 
mass fraction contours. Figure 8.19a shows the x-t plot for water mass fraction for Mach 
5.11 using shock capturing method whereas Fig. 819b shows the same using shock- fitting 
method. It is clear from Fig. 8.19b that almost all the intricate details of the flow field has 
been simulated very clearly whereas in Fig. 8.19a all the flow features are smeared. The 
cellular structure of the detonation wave phenomenon is clearly evident in Fig. 819b. It 
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must be emphasized that the number of grid points used in shock capturing were almost 
three times more than that in shock- fitting, yet the resolution of the flow field was much 
better in shock-fitting method. 

8.4 Ruegg and Dorsey’s Simulation 

In the present study Ruegg and Dorsey’s [6] three cases of ballistic data were 
simulated with the free-stream conditions given in Table 2.2. The grid used in the 
calculations was with 101 grid points in the x direction and 91 points in the y direction. 
This grid was selected based on our previous work [25], in which an extensive grid- 
refinement study was done to determine the grid most suited to this type of study and 
to give a grid-independent solution. 

Figure 8.20a shows the contour plot for pressure for the steady case. The detailed 
mechanism of detonation-wave development depends primarily on the action of pressure 
waves that are generated in the course of combustion. This point has been given special 
attention by earlier researchers [12], Under the current free-stream conditions, nearly 
no pressure waves are generated between the bow shock and the blunt body because of 
the weak reaction wave. As we move to Fig. 8.20b, which corresponds to the regular- 
regime free-stream conditions with p = 0.25 atm, we can see the pressure waves that are 
generated by the combustion between the bow shock and the body. Figure 8.20c shows 
the pressure contours for the large -disturbance-regime case, where p = 0.5 atm. The 
compression waves for this case are much more intense and the bow shock is distorted 
significantly. Inspection of Figs. 8.20a-c shows that the combustion moves the shock 
wave away from the body; this effect increases as the pressure increases. 

Figures 8.21a-c show the density contours for the three cases. The flow has an 
outer bow shock, followed by an induction zone and a reaction front. The induction zone 
(i.e., the separation between the bow shock and the reaction front) is maximum for the 
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steady-regime case when p = 0.1 atm. Both the bow shock and the reaction fronts are 
very smooth in structure. If the pressure rises to 0.25 atm, then periodic instabilities arise 
near the stagnation line and spread downstream, as shown in Fig. 8.21b. Furthermore, 
the induction zone is considerably shorter than the previous case. A jump occurs in the 
density after the bow shock, followed by a drop in density after the reaction front due to 
the increase in temperature. Figure 8.21c shows the density contours when the pressure 
is increased further to 0.5 atm. The induction zone is further decreased as the combustion 
becomes more intense which causes the bow shock to distort. 

Figures 8.22a-c show the computed shadowgraphs for the three cases. In Fig. 8.22a 
both the bow shock and the reaction front are separated considerably, and both fronts 
are smooth. The combustion front is quite weak, as has been observed experimentally. 
Figure 8.22b clearly shows the regular pattern of the instabilities of the reaction front. 
The instabilities are almost periodic in nature. The computed shadowgraph for the large- 
disturbance-regime case is shown in Fig. 8.22c. The distorted bow shock and an irregular, 
unsteady reaction front are clearly evident. A side-by-side comparison of Figs. 8.22b and 
8.22c clearly shows that almost four to five periods of oscillation of the regular regime 
are equivalent to one period of oscillation for the large-disturbance regime. This result 
was reported in Ruegg and Dorsey’s work also [6]. Thus, the large-disturbance regime 
consists of long period oscillations as compared with the regular regime. 

Figures 8.23a-c show a schematic of an oblique detonation wave-engine and three 
possible time history plots of pressure depending upon the various free stream conditions. 
Since this phenomenon has been observed in detonations in tubes as well as in external 
flows, therefore, considering that the same free stream conditions can also be encountered 
in the combustors of an oblique detonation wave engine, and, therefore, similar flow 
features can be expected in the combustors of an oblique detonation wave engine. Figure 
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8.23a shows the x-t plot for pressure when the free stream pressure is 0.1 atm. It shows a 
very smooth bow shock and reaction front. Figure 8.23b shows the x-t plot for pressure 
when the free stream pressure is 0.25 atm. Now the same smooth reaction front has 
turned into a highly periodic reaction front. Compression waves generated between the 
reaction front and the blunt projectile starts interacting with each other and this interaction 
of the compression waves sustains the instability. Moving to Fig. 8.23c where the free 
stream pressure is now 0.5 atm, we notice that the regular periodicity of the reaction 
front disappears and its place is taken over by bumps in the reaction front and the bow 
shock. Both the reaction front and the bow shock are highly distorted. This is the large- 
disturbance regime. 

Figure 8.24a-c show the time history plots for the three cases for water mass fraction 
along the stagnation streamline. Figure 8.24a shows a smooth and steady bow shock and 
reaction front. The maximum water production after the reaction front shows that the 
reaction is nearly completed after the reaction front. As the pressure is increased to 
0.25 atm and Mach number is increased to 4.9, a highly periodic unsteady reaction 
front develops (Fig. 8.24b). The induction zone has reduced considerably. Although the 
reaction front is unsteady, the bow shock is quite smooth. This effect has been observed 
experimentally also. Figure 8.24c shows the large-disturbance-regime case, in which 
the combustion becomes so intense that the bow shock is completely distorted. Highly 
pronounced combustion surges cause a completely different profile of the reaction front 
than for the regular regime. The induction zone is much shorter than for the previous 
two cases. 

In order to validate the present numerical results, we compare the wave-detachment 
distance (both for the shock wave and the combustion wave) for the three cases with the 
experimental results of Ruegg and Dorsey. Figure 8.25 shows the plot of wave detachment 
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versus Mach number for p = 0.1 and 0.5 atm. For both cases, the numerical shock stand- 
off distance and the reaction-front stand-off distance closely match the experimental 
data. Figure 8.26 shows the same plot for p = 0.25 atm. In this case, the numerical 
reaction front closely matches the experimental value, but the shock stand-off distance is 
underpredicted by about 8-9 percent. This discrepancy has not been resolved at present. 

Figure 8.27a shows the time history plot for pressure, and Fig. 8.27c shows the time 
history plot for density for the regular-regime case (i.e, when p = 0.25 atm). This pattern 
of instabilities (which has also been observed in Lehr’s work [11] has been explained in 
the past by many researchers with McVey and Toong’s one-dimensional wave-interaction 
model. A brief description of this model and a discussion of regular-regime results in 
regard to the model has already been explained in the earlier sections. It need only to 
be emphasized here that the results presented in Fig. 8.27a-c are very similar to Lehr’s 
Mach 5.11 case. 

Now we concentrate our attention on the large-disturbance regime (i.e the free- 
stream pressure is increased to 0.5 atm). Figure 8.28a shows the time history plot along 
the stagnation streamline for pressure. Figure 8.28b shows the one-dimensional wave- 
interaction model due to Matsuo et al. [24], and Fig. 8.28c shows the time history plot 
for density. To understand the instabilities for the large-disturbance regime, we must 
consider Figs. 8.28a-c together. The concepts developed in Ref. [12] and used in the 
model by Matsuo et al. [24] are shown in these figures. As shown in Fig. 8.28b, the 
cycle begins when an “explosion within an explosion” occurs on the reaction front; then, 
the reaction shocks propagate upstream and downstream. The forward shock which has 
been referred to as “superdetonation”, moves into the unbumed gases. The backward 
shock, referred to as “retonation”, moves into the burned gases. At a later time, the 
detonation wave overtakes and penetrates the bow shock. Then, the bow shock and 
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the detonation wave create a triple point and generate a reflected shock and a contact 
discontinuity. The bow shock is accelerated by this penetration, and the gas behind 
the bow shock is further compressed. At a later time, the bow shock is decelerated, 
and the bow shock and reaction front become separated. The retonation wave and the 
reflected shock reach the body surface, and the reflected compression wave goes to the 
bow shock. The reflected compression wave interacts with the bow shock, and another 
contact discontinuity is generated. Finally, the contact discontinuity reaches the original 
reaction front, the explosion within an explosion occurs on reaction front, and the cycle 
of events is repeated. 

By comparing the model (Fig. 8.28b) with the numerical results (Figs. 8.28a and 
8.28c), we see that the model developed by Matsuo et al. [24] agrees well with the current 
numerical results. As the pressure is increased to 0.5 atm, the strong exothermicity on 
the reaction front generates strong reaction shocks, which move toward the bow shock 
and the body. The exothermic reaction accelerates the reaction shock. As mentioned 
earlier, this phenomenon is called the explosion within an explosion and has also been 
observed in detonations in tubes [12]. Both Figs. 8.28a and 8.28c clearly show this 
phenomenon. At this point explosion, two strong shocks are produced one that moves 
into the unbumed medium, which we refer to as “superdetonation,” and the other that 
moves into the burned gases and is known as “retonation.” Both superdetonation and 
retonation are seen in the numerical results of Fig. 8.28a and Fig. 8.28c. 

So that we could better examine the series of wave interactions and the many point 
explosions and penetrations, the large-disturbance case was run longer. Figures 8.29a 
and 8.29b show the time history plots for pressure and density for an extended period 
of time. Two penetrations are evident (i.e, transitions from deflagration to detonation 
and vice versa). 
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Figures 8.30a and 8.30b show enlarged views of the history plots for pressure 
and density. As shown in Fig. 8.31a, the superdetonation causes the reaction front to 
accelerate, which in turn causes the bulging, or penetration, of the bow shock. This 
case is a typical example of deflagration-to-detonation (DDT) transition and also has 
been observed in detonations in tubes. Because of the accelerating reaction front and the 
eventual merging of the bow shock with the reaction front, a self "Sustained detonation 
wave exists in front of the projectile body. The point of intersection of the bow shock 
and reaction front is termed as “triple point”. A reflected shock is necessary to satisfy the 
continuity of the pressure and flow direction across the point of intersection of the bow 
shock and reaction front. This reflected shock is clearly shown in the numerical results. 
Thus, the triple point is the point at which the bow shock, the reaction front and the 
reflected shock meet. Of the three phenomena, the reflected shock is the weakest. Next, 
the bow shock wave begins to decelerate. Thus, the merged bow shock and reaction front 
begin to separate, as shown in Fig. 8.30b. This separation stage is the transition from 
detonation to deflagration; thus, it is known as the ordinary shock-induced combustion. 
The reflected shock (after reflecting from the body) reaches the bow shock and generates 
a contact discontinuity (as shown in Fig. 8.30b), and the entire process is repeated. 
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Figure 8.20 Pressure contours for (a) Stable case, (b) Unstable regular regime, and (c) Unstable large-disturbance regime. 
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Figure 8.21 Density contours for (a) Stable case, (b) Unstable regular regime, and (c) Unstable large-disturbance regime. 




Figure 8.22 Computed shadowgraphs for (a) Stable case, (b) Unstable regular regime, and (c) Unstable 

large-disturbance regime. 









TIME HISTORY PLOTS OF WATER MASS FRACTION 



ST^ADYJ-jEGIMEjP_jJDJ_atm)^^^^REGULAFnjNSTEA^ UNSTEADY REGIME (P = 0.5 atm) 


Figure 8.24 Time history plots of water mass fraction along stagnation line for (a) Stable case, (b) Unstable regular 

regime, and (c) Unstable large-disturbance regime. 
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Figure 8.25 Comparison of wave detachment distance from experiments of Ruegg and Dorsey, and 
the current numerical work for pressure of 0.1 atm and 0.5 atm. 
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Figure 8.26 Comparison of wave detachment distance from experiments of Ruegg and Dorsey, and the 

current numerical work for pressure of 0.25 atm. 




Figure 8.27 Side-by-side comparison of time history plots for regular regime and the wave interaction model: (a) Pressure 
contours from simulation, (b) McVey and Toong wave interaction model, and (c) Density contours from simulation. 
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Figure 8.28 Side-by-side comparison of time history plots for large-disturbance regime and the wave interaction model: (a) 
Pressure contours from simulation, (b) Matsuo et al. wave interaction model, and (c) Density contours from simulation. 
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Figure 8.30 Enlarged view of x-t plots for large-disturbance regime showing DDT; (a) Pressure contours from 

simulation and (b) Density contours from simulation. 









Chapter 9 


CONCLUSIONS 

A numerical study is presented to demonstrate the capabilities of modem CFD tech- 
niques for complex reacting flow-field predictions. Unsteady and steady ballistic-range 
experiments are successfully simulated. The calculations have been carried out for Mach 
5.11 and 6.46. The Mach 5.11 case was found to be unsteady with periodic oscillations. 
The frequency of oscillations was calculated and was found to be in good agreement 
with the experimentally observed frequency. The Mach 6.46 case was found to be 
macroscopically stable, thus supporting the existing view that it is possible to stabilize 
the shock-induced combustion phenomena with sufficient level of overdrive. All the ex- 
perimentally observed flow features are captured, and the frequency of the combustion 
instabilities is found to be in good agreement with the experimentally observed frequency. 

Study also included simulations of three cases of Ruegg and Dorsey’s ballistic-range 
experiments. Results show that an increase in the free-stream pressure can drive a sta- 
ble reaction front to a regular, periodic unstable regime and again to a large-disturbance 
regime, as observed experimentally. The flow features of the regular regime have been 
successfully described with the one-dimensional wave-interaction model of McVey and 
Toong. The flow features of the large-disturbance regime have been described by the 
one-dimensional model of Matsuo et al. The shock stand-off and combustion stand-off 
distances agree well with the experimental results of Ruegg and Dorsey. 

The current study used both shock-capturing and shock-fitting techniques to simu- 
late the shock-induced combustion past blunt projectiles. In such problems which involve 


158 



instabilities, the shock-fitting technique gives much better resolution of the flow features 
than the shock-capturing technique. The observed flow features have been successfully 
correlated with the one-dimensional wave-interaction models, and the frequency of oscil- 
lations has been matched with the experimental data as well as with earlier investigations. 
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APPENDIX A 


SOLUTION OF GOVERNING EQUATIONS ALONG STAGNATION 
LINE FOR AXISYMMETRIC FLOWS USING L’HOPITAL’S RULE 


L’hopital rule says that if the function fix) and <f>(x) vanish at the point x = a, that 
i sfia) = <j>(a) = 0, then if the ratio has a limit a s x —> a there also exists a limit 
lim^fj and = lim^fj . Since for axisymmmetric problems we have a 

coordinate singularity at the symmetry axis, y = 0, in the source term, therefore, to 
remove the singularity we use L’Hopital’s rule and make the following replacements in 
the equations: 


V 

dv 

55 

= =1 

y 

§? 

1 

Vx_ 

dv x 
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_ ^iy 

y 

if 

1 

% 
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~5y 

_ _ 


^ it 1 


n fv v yy + v y- v y\ 

~ 2/X V 2 y ) ~ ^ Vyy 



k OT = kTy = = kT 

y dy y | J 


With these changes the governing equations at the axis of symmetry reduces to: 
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u = 


pU^ + p + Txx 
pUV + Tyx 

_(pE + P 4 ■ T xx )u + T yx v+ q x 
pv 

pVU + Tyx 

pv 2 + P + Tyy 

_(pE “I - p “I - 4~ Tyx u4"C[j/ 

P v y 

PUVy-p[V X y + Uyy] 
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pEVy + pu y - ~p[2v 2 y - Vy - Vy U* ] ~ ^ [ V X y U + U Uyy ] + liTj 


(See Below) 
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Also the following terms to be replaced with their final values given on the R.H.S. 


2 / du dv v 

xx S^\dx dy y 


2 / du dv dv 

3 ^ \ dx dy dy 


= --[1(2— — 2 —' 
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Hence, 


and 


2 / dv v du\ 

t « = _ 3 l ‘Ydj~y-te) 

2 / dv dv du\ 

3 \ dy dy dx) 


2 f dv du\ 

Tyy = ~V\di~fo) 


2 / v dv du\ 
T 6 e = -^[ 2 y-fy-fc) 

2 / dv dv du\ 
3^\ dy dy dx) 

2 f dv du\ 

T ee = - 3 fi [dy~dx) 


Though the governing equations changes at the axis of symmetry as shown above, 
but we would still like to use the same set of governing equations everywhere (including 
at the axis of symmetry) with only a slight modifications at the axis of symmetry to 
take care of the above changes. The only term needs to be taken care of at the axis of 
symmetry is the term ^ 

As y — * 0 and v — ► 0 at the axis of symmetry we get ^ = indeterminate. 

Therefore, using L’Hopital’s rule we have 

dv a 

- - N. - — 

y ~ % " a y 

_ V(2) — U(l) _ U(2) - 0 _ V(2) 

~ y{ 2 ) - y{ i) _ y( 2 ) - o y( 2 ) 

(Values inside the brackets denote grid points along the body) 

Hence, 

V t>(2) 

y 2/(2) 
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This is because both the normal velocity and value of y is zero at the stagnation line. 
Therefore, at the axis of symmetry simply replace J terms by gf} terms. Also if absolute 
values of y(i,i) is > 10 -10 (machine zero) then use the general expression for | i.e., at 
the stagnation line it will be 

But if the absolute value of y(i, 1 ) is < 10~ 10 which means zero, then L’Hopital’s 
rule is used for ^ i.e., ^ i.e., the value at the second grid point is used at the 

stagnation line. 
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APPENDIX B 


CALCULATION OF THE CHEMICAL JACOBIAN MATRIX 

We shall be discussing the calculations of the chemical jacobian matric for the species 
0 2 , H 2 0, H 2 , OH 

Applying the law of mass action to the global model 2H 2 + 0 2 2H 2 0 we have 


Co 2 = — kf 2 C^C 0i + fcftiCoH 

(Bl) 

Ch 2 o = 2\k h C 2 ou C E2 - Ch 20 ] 
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Ch 2 = Co* ~ 2 ^ 2 ° 
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Con = - (2C02 + C'hjo) 

(B4) 
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2 — > 0 2 


OH 


4 -> H 2 0 


Derivation of u>,- terms in terms of U Variables 

Note: This is required before we carry out differentiation i.e. ^ r e ) 

<7U(Ui ,U2, Us) * 


Now, 
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Therefore, 
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if density is given in and molecular weight in gra m-mole 


Therefore, 
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From Eq. (Bl) we have 
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From Eq. (B3) we have 



Similarly, from Eq. (B2) we have, 
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